Marie Chavent
email: marie.chavent@math.u-bordeaux.fr

Vanessa Kuentz-Simonet
email: vanessa.kuentz-simonet@irstea.fr

Amaury Labenne
email: amaury.labenne@yahoo.fr

Jérôme Saracco
email: jerome.saracco@ensc.fr

ClustGeo: an R package for hierarchical clustering with spatial constraints

Keywords: Ward-like hierarchical clustering, Soft contiguity constraints, Pseudo-inertia, Non-Euclidean dissimilarities, Geographical distances

In this paper, we propose a Ward-like hierarchical clustering algorithm including spatial/geographical constraints. Two dissimilarity matrices D 0 and D 1 are inputted, along with a mixing parameter α ∈ [0, 1]. The dissimilarities can be non-Euclidean and the weights of the observations can be non-uniform. The first matrix gives the dissimilarities in the "feature space" and the second matrix gives the dissimilarities in the "constraint space". The criterion minimized at each stage is a convex combination of the homogeneity criterion calculated with D 0 and the homogeneity criterion calculated with D 1 . The idea is then to determine a value of α which increases the spatial contiguity without deteriorating too much the quality of the solution based on the variables of interest i.e. those of the feature space. This procedure is illustrated on a real dataset using the R package ClustGeo.

Introduction

The difficulty of clustering a set of n objects into k disjoint clusters is one that is well known among researchers. Many methods have been proposed either to find the best partition according to a dissimilarity-based homogeneity criterion, or to fit a mixture model of multivariate distribution function. However, in some clustering problems, it is relevant to impose constraints on the set of allowable solutions. In the literature, a variety of different solutions have been suggested and applied in a number of fields, including earth science, image processing, social science, and -more recently -genetics. The most common type of constraints are contiguity constraints (in space or in time). Such restrictions occur when the objects in a cluster are required not only to be similar to one other, but also to comprise a contiguous set of objects. But what is a contiguous set of objects?

Consider first that the contiguity between each pair of objects is given by a matrix C = (c ij ) n×n , where c ij = 1 if the ith and the jth objects are regarded as contiguous, and 0 if they are not. A cluster C is then considered to be contiguous if there is a path between every pair of objects in C (the subgraph is connected). Several classical clustering algorithms have been modified to take this type of constraint into account (see e.g. [START_REF] Murtagh | Multidimensional clustering algorithms[END_REF], [START_REF] Legendre | [END_REF], [4]). Surveys of some of these methods can be found in [START_REF] Gordon | A survey of constrained classification[END_REF] and [START_REF] Murtagh | A Survey of Algorithms for Contiguity-constrained Clustering and Related Problems[END_REF]. For instance, the standard hierarchical procedure based on Lance and Williams formula [10] can be constrained by merging only contiguous clusters at each stage. But what defines "contiguous" clusters? Usually, two clusters are regarded as contiguous if there are two objects, one from each cluster, which are linked in the contiguity matrix. But this can lead to reversals (i.e. inversions, upward branching in the tree) in the hierarchical classification. It was proven that only the complete link algorithm is guaranteed to produce no reversals when relational constraints are introduced in the ordinary hierarchical clustering procedure [START_REF] Ferligoj | Clustering with relational constraint[END_REF]. Recent implementation of strict constrained clustering procedures are available in the R package const.clust [START_REF] Legendre | const.clust: Space-and Time-Constrained Clustering Package[END_REF] and in the python library clusterpy1 . Hierarchical clustering of SNPs (Single Nucleotide Polymorphism) with strict adjacency constraint is also proposed in [7] and implemented in the R package BALD 2 . The recent R package Xplortext [START_REF] Bécue-Bertaut | Xplortext: Statistical Analysis of Textual Data R package[END_REF] implements also chronogically constrained agglomerative hierarchical clustering for the analysis of textual data.

The previous procedures which impose strict contiguity may separate objects which are very similar into different clusters, if they are spatially apart. Other non-strict constrained procedures have then been developed, including those referred to as soft contiguity or spatial constraints. For example, Oliver et al. [16] and Bourgault et al. [START_REF] Bourgault | The Multivariate (co) Variogram as a Spatial Weighting Function in Classification Methods[END_REF] suggest running clustering algorithms on a modified dissimilarity matrix . This dissimilarity matrix is a combination of the matrix of geographical distances and the dissimilarity matrix computed from non geographical variables. According to the weights given to the geo-graphical dissimilarities in this combination, the solution will have more or less spatially contiguous clusters. However, this approach raises the problem of defining weight in an objective manner.

In image processing, there are many approaches for image segmentation including for instance usage of convolution and wavelet transforms. In this field non-strict spatially constrained clustering methods have been also developed. Objects are pixels and the most common choices for the neighborhood graph are the four and eight neighbors graphs. A contiguity matrix C is used (and not a geographical dissimilarity matrix as previously) but the clusters are not strictly contiguous, as a cluster of pixels does not necessarily represent a single region on the image. Ambroise et al. [START_REF] Ambroise | Clustering of spatial data by the EM algorithm[END_REF] [2] suggest a clustering algorithm for Markov random fields based on an EM (Expectation-Maximization) algorithm. This algorithm maximizes a penalized likelihood criterion and the regularization parameter gives more or less weight to the spatial homogeneity term (the penalty term). Recent implementations of spatially-located data clustering algorithms are available in SpaCEM3 3 , dedicated to Spatial Clustering with EM and Markov Models. This software uses the model proposed in [START_REF] Vignes | Gene Clustering via Integrated Markov Models Combining Individual and Pairwise Features[END_REF] for gene clustering via integrated Markov models. In a similar vein, Miele et al. [13] proposed a model-based spatially constrained method for the clustering of ecological networks. This method embeds geographical information within a EM regularization framework by adding some constraints to the maximum likelihood estimation of parameters. The associated R package is available at http://lbbe.univ-lyon1.fr/geoclust. Note that all these methods are partitioning methods and that the constraints are neighborhood constraints.

In this paper, we propose a hierarchical clustering (and not partitioning) method including spatial constraints (not necessarily neighborhood constraints). This Ward-like algorithm uses two dissimilarity matrices D 0 and D 1 and a mixing parameter α ∈ [0, 1]. The dissimilarities are not necessarily Euclidean (or non Euclidean) distances and the weights of the observations can be non-uniform. The first matrix gives the dissimilarities in the 'feature space' (socio-economic variables or grey levels for instance). The second matrix gives the dissimilarities in the 'constraint space'. For instance, D 1 can be a matrix of geographical distances or a matrix built from the contiguity matrix C. The mixing parameter α sets the importance of the constraint in the clustering procedure. The criterion minimized at each stage is a convex combination of the homogeneity criterion calculated with D 0 and the homogeneity criterion calculated with D 1 . The parameter α (the weight of this convex combination) controls the weight of the constraint in the quality of the solutions. When α increases, the homogeneity calculated with D 0 decreases whereas the homogeneity calculated with D 1 increases. The idea is to determine a value of α which increases the spatial-contiguity without deteriorating too much the quality of the solution on the variables of interest. The R package ClustGeo [START_REF] Chavent | ClustGeo: Hierarchical Clustering with Spatial Constraints[END_REF] implements this constrained hierarchical clustering algorithm and a procedure for the choice of α.

The paper is organized as follows. After a short introduction (this section), Section 2 presents the criterion optimized when the Lance-Williams parameters [10] as used in Ward's minimum variance method but dissimilarities are not necessarily Euclidean (or non-Euclidean) distances. We also show how to implement this procedure with the package ClustGeo (or the R function hclust) when non-uniform weights are used. In Section 3 we present the constrained hierarchical clustering algorithm which optimizes a convex combination of this criterion calculated with two dissimilarity matrices. Then the procedure for the choice of the mixing parameter is presented as well as a description of the functions implemented in the package ClustGeo. In Section 4 we illustrate the proposed hierarchical clustering process with geographical constraints using the package ClustGeo before a brief discussion given in Section 5.

Throughout the paper, a real dataset is used for illustration and reproducibility purposes. This dataset contains 303 French municipalities described based on four socio-economic variables. The matrix D 0 will contain the socioeconomic distances between municipalities and the matrix D 1 will contain the geographical distances. The results will be easy to visualize on a map.

Ward-like hierarchical clustering with dissimilarities and non-uniform weights

Let us consider a set of n observations. Let w i be the weight of the ith observation for i = 1, . . . , n. Let D = [d ij ] be a n × n dissimilarity matrix associated with the n observations, where d ij is the dissimilarity measure between observations i and j. Let us recall that the considered dissimilarity matrix D is not necessarily a matrix of Euclidean (or non Euclidean) distances. When D is not a matrix of Euclidean distances, the usual inertia criterion (also referred to as variance criterion) used in Ward hierarchical clustering approach [START_REF] Ward | Hierarchical Grouping to Optimize an Objective Function[END_REF] is meaningless and the Ward algorithm implemented with the Lance and Williams [10] formula has to be re-interpreted. The Ward method has already been generalized to use with non Euclidean distances, see e.g. [START_REF] Strauss | ward?s method for use with manhattan distances[END_REF] for l 1 norm or Manhattan distances. In this section the more general case of dissimilarities is studied. We first present the homogeneity criterion which is optimized in that case and the underlying aggregation measure which leads to a Wardlike hierarchical clustering process. We then provide an illustration using the package ClustGeo and the well-known R function hclust.

The Ward-like method

Pseudo-inertia. Let us consider a partition

P K = (C 1 , . . . , C K ) in K clusters.
The pseudo-inertia of a cluster C k generalizes the inertia to the case of dissim-Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1 ilarity data (Euclidean or not) in the following way :

I(C k ) = i∈C k j∈C k w i w j 2µ k d 2 ij (1)
where µ k = i∈C k w i is the weight of C k . The smaller the pseudo-inertia I(C k ) is, the more homogenous are the observations belonging to the cluster C k . The pseudo within-cluster inertia of the partition P K is therefore:

W (P K ) = K k=1 I(C k ).
The smaller this pseudo within-inertia W (P K ) is, the more homogenous is the partition in K clusters.

Spirit of the Ward hierarchical clustering. To obtain a new partition P K in K clusters from a given partition P K+1 in K + 1 clusters, the idea is to aggregate the two clusters A and B of P K+1 such that the new partition has minimum within-cluster inertia (heterogeneity, variance), that is:

arg min A,B∈P K+1 W (P K ), (2)
where P K = P K+1 \{A, B} ∪ {A ∪ B} and

W (P K ) = W (P K+1 ) -I(A) -I(B) + I(A ∪ B).
Since W (P K+1 ) is fixed for a given partition P K+1 , the optimization problem (2) is equivalent to:

min A,B∈P K+1 I(A ∪ B) -I(A) -I(B). (3)
The optimization problem is therefore achieved by defining

δ(A, B) := I(A ∪ B) -I(A) -I(B)
as the aggregation measure between two clusters which is minimized at each step of the hierarchical clustering algorithm. Note that δ(A, B) = W (P K ) -W (P K+1 ) can be seen as the increase of within-cluster inertia (loss of homogeneity).

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

Ward-like hierarchical clustering process for non-Euclidean dissimilarities.

The interpretation of the Ward hierarchical clustering process in the case of dissimilarity data is the following:

-Step K = n: initialization. The initial partition P n in n clusters (i.e. each cluster only contains an observation) is unique. -Step K = n -1, . . . , 2: obtaining the partition in K clusters from the partition in K + 1 clusters. At each step K, the algorithm aggregates the two clusters A and B of P K+1 according to the optimization problem (3) such that the increase of the pseudo within-cluster inertia is minimum for the selected partition over the other ones in K clusters. -Step K = 1: stop. The partition P 1 in one cluster (containing the n observations) is obtained.

The hierarchically-nested set of such partitions {P n , . . . , P K , . . . , P 1 } is represented graphically by a tree (also called dendrogram) where the height of a cluster

C = A ∪ B is h(C) := δ(A, B).
In practice, the aggregation measures between the new cluster A ∪ B and any cluster D of P K+1 are calculated at each step thanks to the well-known Lance and Williams [10] equation:

δ(A ∪ B, D) = µ A + µ D µ A + µ B + µ D δ(A, D) + µ B + µ D µ A + µ B + µ D δ(B, D) - µ D µ A + µ B + µ D δ(A, B). (4)
In the first step the partition is P n and the aggregation measures between the singletons are calculated with

δ ij := δ({i}, {j}) = w i w j w i + w j d 2 ij ,
and stored in the n × n matrix ∆ = [δ ij ]. For each subsequent step K, the Lance and Williams formula (4) is used to build the corresponding K × K aggregation matrix. The hierarchical clustering process described above is thus suited for non-Euclidean dissimilarities and then for non-numerical data. In this case, it optimises the pseudo within-cluster inertia criterion [START_REF] Bécue-Bertaut | Xplortext: Statistical Analysis of Textual Data R package[END_REF].

Case when the dissimilarities are Euclidean distances. When the dissimilarities are Euclidean distances calculated from a numerical data matrix X of dimension n × p for instance, the pseudo-inertia of a cluster C k defined in ( 1) is now equal to the inertia of the observations in C k :

I(C k ) = i∈C k w i d 2 (x i , g k )
Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

where x i ∈ p is the ith row of X associated with the ith observation, and

g k = 1 µ k i∈C k w i x i ∈ R p is the center of gravity of C k .
The aggregation measure δ(A, B) between two clusters is written then as:

δ(A, B) = µ A µ B µ A + µ B d 2 (g A , g B ).

Illustration using the package ClustGeo

Let us examine how to properly implement this procedure with R. The dataset is made up of n = 303 French municipalities described based on p = 4 quantitative variables and is available in the package ClustGeo. A more complete description of the data is provided in Section 4.1.

> library(ClustGeo) > data(estuary) > names(estuary) [START_REF] Ambroise | Clustering of spatial data by the EM algorithm[END_REF] "dat" "D.geo" "map"

To carry out Ward hierarchical clustering, the user can use the function hclustgeo implemented in the package ClustGeo taking the dissimilarity matrix D (which is an object of class dist, i.e. an object obtained with the function dist or a dissimilarity matrix transformed in an object of class dist with the function as.dist) and the weights w = (w 1 , . . . , w n ) of observations as arguments.

> D <-dist(estuary$dat) > n <-nrow(estuary$dat) > tree <-hclustgeo(D,wt=rep(1/n,n))
Remarks.

-The function hclustgeo is a wrapper of the usual function hclust with the following arguments:

method = "ward.D", d = ∆, members = w.

For instance, when the observations are all weighted by 1/n , the argument d must be the matrix ∆ = D 2 2n and not the dissimilarity matrix D:

> tree <-hclust(D^2/(2*n),method="ward.D")
-As mentioned before, the user can check that the sum of the heights in the dendrogram is equal to the total pseudo-inertia of the dataset:

> inertdiss(D,wt=rep(1/n,n)) # the pseudo-inertia of the data [1] 1232.769 > sum(tree$height) [1] 1232.769
-When the weights are not uniform, the calculation of the matrix ∆ takes a few lines of code and the use of the function hclustgeo is clearly more convenient than hclust:

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

> w <-estuary$map@data$POPULATION # non-uniform weights > tree <-hclustgeo(D,wt=w) > sum(tree$height) [START_REF] Ambroise | Clustering of spatial data by the EM algorithm[END_REF] 1907989

versus > Delta <-D > for (i in 1:(n-1)) { for (j in (i+1):n) { Delta[n*(i-1) -i*(i-1)/2 + j-i] <- Delta[n*(i-1) -i*(i-1)/2 + j-i]^2*w[i]*w[j]/(w[i]+w[j]
)}} > tree <-hclust(Delta,method="ward.D",members=w) > sum(tree$height) [START_REF] Ambroise | Clustering of spatial data by the EM algorithm[END_REF] 1907989 [START_REF] Bécue-Bertaut | Xplortext: Statistical Analysis of Textual Data R package[END_REF]

Ward-like hierarchical clustering with two dissimilarity matrices

Let us consider again a set of n observations, and let w i be the weight of the ith observation for i = 1, . . . , n. Let us now consider that two n × n dissimilarity matrices D 0 = [d 0,ij ] and D 1 = [d 1,ij ] are provided. For instance, let us assume that the n observations are municipalities, D 0 can be based on a numerical data matrix of p 0 quantitative variables measured on the n observations and D 1 can be a matrix containing the geographical distances between the n observations. In this section, a Ward-like hierarchical clustering algorithm is proposed. A mixing parameter α ∈ [0, 1] allows the user to set the importance of each dissimilarity matrix in the clustering procedure. More specifically, if D 1 gives the dissimilarities in the constraint space, the mixing parameter α sets the importance of the constraint in the clustering procedure and controls the weight of the constraint in the quality of the solutions.

Hierarchical clustering algorithm with two dissimilarity matrices

For a given value of α ∈ [0, 1], the algorithm works as follows. Note that the partition in K clusters will be hereafter indexed by α as follows: P α K .

Definitions. The mixed pseudo inertia of the cluster C α k (called mixed inertia hereafter for sake of simplicity) is defined as

I α (C α k ) = (1 -α) i∈C α k j∈C α k w i w j 2µ α k d 2 0,ij + α i∈C α k j∈C α k w i w j 2µ α k d 2 1,ij , (5)
where

µ α k = i∈C α k w i is the weight of C α k and d 0,ij (resp. d 1,ij
) is the normalized dissimilarity between observations i and j in D 0 (resp. D 1 ).

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

The mixed pseudo within-cluster inertia (called mixed within-cluster inertia hereafter for sake of simplicity) of a partition P α K = (C α 1 , . . . , C α K ) is the sum of the mixed inertia of its clusters:

W α (P α K ) = K k=1 I α (C α k ). ( 6
)
Spirit of the Ward-like hierarchical clustering. As previously, in order to obtain a new partition P α K in K clusters from a given partition P α K+1 in K +1 clusters, the idea is to aggregate the two clusters A and B of P K+1 such that the new partition has minimum mixed within-cluster inertia. The optimization problem can now be expressed as follows:

arg min A,B∈P α K+1 I α (A ∪ B) -I α (A) -I α (B). ( 7
)
Ward-like hierarchical clustering process.

-Step K = n: initialization. The dissimilarities can be re-scaled between 0 and 1 to obtain the same order of magnitude: for instance,

D 0 ← D 0 max(D 0 )
and

D 1 ← D 1 max(D 1 )
.

Note that this normalization step can also be done in a different way.

The initial partition P α n =: P n in n clusters (i.e. each cluster only contains an observation) is unique and thus does not depend on α.

-Step K = n -1, . . . , 2: obtaining the partition in K clusters from the partition in K + 1 clusters. At each step K, the algorithm aggregates the two clusters A and B of P α K+1 according to the optimization problem (7) such that the increase of the mixed within-cluster inertia is minimum for the selected partition over the other ones in K clusters. More precisely, at step K, the algorithm aggregates the two clusters A and B such that the corresponding aggregation measure

δ α (A, B) := W α (P α K+1 ) -W α (P α K ) = I α (A ∪ B) -I α (A) -I α (B)
is minimum.

-Step K = 1: stop. The partition P α 1 =: P 1 in one cluster is obtained. Note that this partition is unique and thus does not depend on α.

In the dendrogram of the corresponding hierarchy, the height of a cluster A∪B is given by δ α (A, B).

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

In practice, the Lance and Williams equation ( 4) remains true in this context (where δ must be replaced by δ α ). The aggregation measure between two singletons are written now:

δ α ({i}, {j}) = (1 -α) w i w j w i + w j d 2 0,ij + α w i w j w i + w j d 2 1,ij .
The Lance and Williams equation is then applied to the matrix

∆ α = (1 -α)∆ 0 + α∆ 1 .
where ∆ 0 (resp.

∆ 1 ) is the n × n matrix of the values δ 0,ij = wiwj wi+wj d 2 0,ij (resp. δ 1,ij = wiwj wi+wj d 2 1,ij ).
Remarks.

-The proposed procedure is different from applying directly the Ward algorithm to the "dissimilarity" matrix obtained via the convex combination

D α = (1-α)D 0 +αD 1 .
The main benefit of the proposed procedure is that the mixing parameter α clearly controls the part of pseudo inertia due to D 0 and D 1 in (5). This is not the case when applying directly the Ward algorithm to D α since it is based on a unique pseudo inertia.

-When α = 0 (resp. α = 1), the hierarchical clustering is only based on the dissimilarity matrix D 0 (resp. D 1 ). A procedure to determine a suitable value for the mixing parameter α is proposed hereafter, see Section 3.2.

A procedure to determine a suitable value for the mixing parameter α

The key point is the choice of a suitable value for the mixing parameter α ∈ [0, 1]. This parameter logically depends on the number of clusters K and this logical dependence is an issue when it comes to decide an optimal value for the parameter α. In this paper a practical (but not globally optimal) solution to this issue is proposed : conditioning on K and choosing α that best compromises between loss of socio-economic and loss of geographical homogeneity. Of course other solutions than conditioning on K could be explore (conditioning on α or defining a global criterion) but these solutions seems to be more difficult to implement in a sensible procedure.

To illustrate the idea of the proposed procedure, let us assume that the dissimilarity matrix D 1 contains geographical distances between n municipalities, whereas the dissimilarity matrix D 0 contains distances based on a n × p 0 data matrix X 0 of p 0 socio-economic variables measured on these n municipalities. An objective of the user could be to determine a value of α which increases the geographical homogeneity of a partition in K clusters without adversely affecting socio-economic homogeneity. These homogeneities can be measured using the appropriate pseudo within-cluster inertias.

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

Let β ∈ [0, 1]. Let us introduce the notion of proportion of the total mixed (pseudo) inertia explained by the partition P α K in K clusters:

Q β (P α K ) = 1 - W β (P α K ) W β (P 1 ) ∈ [0, 1].
Some comments on this criterion.

-When β = 0, the denominator W 0 (P 1 ) is the total (pseudo) inertia based on the dissimilarity matrix D 0 and the numerator is the (pseudo) withincluster inertia W 0 (P α K ) based on the dissimilarity matrix D 0 , i.e. only from the socio-economic point of view in our illustration. The higher the value of the criterion Q 0 (P α K ), the more homogeneous the partition P α K is from the socio-economic point of view (i.e. each cluster C α k , k = 1, . . . , K has a low inertia I 0 (C α k ) which means that individuals within the cluster are similar). When the considered partition P α K has been obtained with α = 0, the criterion Q 0 (P α K ) is obviously maximal (since the partition P 0 K was obtained by using only the dissimilarity matrix D 0 ), and this criterion will naturally tend to decrease as α increases from 0 to 1.

-Similarly, when β = 1, the denominator W 1 (P 1 ) is the total (pseudo) inertia based on the dissimilarity matrix D 1 and the numerator is the (pseudo) within-cluster inertia W 1 (P α K ) based on the dissimilarity matrix D 1 , i.e. only from a geographical point of view in our illustration. Therefore, the higher the value of the criterion Q 1 (P α K ), the more homogeneous the partition P α K from a geographical point of view. When the considered partition P α K has been obtained with α = 1, the criterion Q 1 (P α K ) is obviously maximal (since the partition P 1 K was obtained by using only the dissimilarity matrix D 1 ), and this criterion will naturally tend to decrease as α decreases from 1 to 0.

-For a value of β ∈]0, 1[, the denominator W β (P 1 ) is a total mixed (pseudo) inertia which can not be easily interpreted in practice, and the numerator W β (P α K ) is the mixed (pseudo) within-cluster inertia. Note that when the considered partition P α K has been obtained with α = β, the criterion Q β (P α K ) is obviously maximal by construction, and it will tend to decrease as α moves away from β.

-Finally, note that this criterion Q β (P α K ) is decreasing in K. Moreover, ∀β ∈ [0, 1], it is easy to see that Q β (P n ) = 1 and Q β (P 1 ) = 0. The more clusters there are in a partition, the more homogeneous these clusters are (i.e. with a low inertia). Therefore this criterion can not be used to select an appropriate number K of clusters.

How to use this criterion to select the mixing parameter α. Let us focus on the above mentioned case where the user is interested in determining a value of α which increases the geographical homogeneity of a partition in K clusters

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

without deteriorating too much the socio-economic homogeneity. For a given number K of clusters (the choice of K is discussed later), the idea is the following:

-Let us consider a given grid of J values for α ∈ [0, 1]:

G = {α 1 = 0, α 1 , . . . , α J = 1}.
For each value α j ∈ G, the corresponding partition P αj K in K clusters is obtained using the proposed Ward-like hierarchical clustering algorithm.

-For the J partitions {P αj K , j = 1, . . . , J}, the criterion Q 0 (P αj K ) is evaluated. The plot of the points {(α j , Q 0 (P αj K )), j = 1, . . . , J} provides a visual way to observe the loss of socio-economic homogeneity of the partition P αj K (from the "pure" socio-economic partition P 0 K ) as α j increases from 0 to 1.

-Similarly, for the J partitions {P αj K , j = 1, . . . , J}, the criterion Q 1 (P αj K ) is evaluated. The plot of the points {(α j , Q 1 (P αj K )), j = 1, . . . , J} provides a visual way to observe the loss of geographical homogeneity of the partition P αj K (from the "pure" geographical partition P 1 K ) as α j decreases from 1 to 0.

-These two plots (superimposed in the same figure) allow the user to choose a suitable value for α ∈ G which is a trade-off between the loss of socioeconomic homogeneity and greater geographical cohesion (when viewed through increasing values of α) .

Case where the two total (pseudo) inertias W 0 (P 1 ) and W 1 (P 1 ) used in Q 0 (P α K ) and Q 1 (P α K ) are very different. Let us consider for instance that the dissimilarity matrix D 1 is a "neighborhood" dissimilarity matrix, constructed from the corresponding adjacency matrix A: that is D 1 = 1 n -A with 1 n,ij = 1 for all (i, j), a ij equal to 1 if observations i and j are neighbors and 0 otherwise, and a ii = 1 by convention. With this kind of local dissimilarity matrix D 1 , the geographical cohesion for few clusters is often small: indeed, W 1 (P 1 ) could be very small and thus the criterion Q 1 (P α K ) takes values generally much smaller than those obtained by the Q 0 (P α K ). Consequently, it is not easy for the user to select easily and properly a suitable value for the mixing parameter α since the two plots are in two very different scales.

One way to remedy this problem is to consider a renormalization of the two plots.

Rather than reasoning in terms of absolute values of the criterion Q 0 (P α K ) (resp. Q 1 (P α K )) which is maximal in α = 0 (resp. α = 1), we will renormalize Q 0 (P α K ) and Q 1 (P α K ) as follows:

Q * 0 (P α K ) = Q 0 (P α K )/Q 0 (P 0 K ) and Q * 1 (P α K ) = Q 1 (P α K )/Q 1 (P 1 K
) and we then reason in terms of proportions of these criteria. Therefore the corresponding plot {(α j , Q * 0 (P αj K )), j = 1, . . . , J} (resp. {(α j , Q * 1 (P αj K )), j = 1, . . . , J}) starts from 100% and decreases as α j increases from 0 to 1 (resp. as α j decreases from 1 to 0).

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

The choice of the number K of clusters. The proposed procedure to select a suitable value for the mixing parameter α works for a given number K of clusters. Thus, it is first necessary to select K. One way of achieving this is to focus on the dendrogram of the hierarchicallynested set of such partitions {P 0 n = P n , . . . , P 0 K , . . . , P 0 1 = P 1 } only based on the dissimilarity matrix D 0 (i.e. for α = 0, that is considering only the socioeconomic point of view in our application). According to the dendrogram, the user can select an appropriate number K of clusters with their favorite rule.

Description of the functions of the package ClustGeo

The previous Ward-like hierarchical clustering procedure is implemented in the function hclustgeo with the following arguments: hclustgeo(D0, D1 = NULL, alpha = 0, scale = TRUE, wt = NULL) where:

-D0 is the dissimilarity matrix D 0 between n observations. It must be an object of class dist, i.e. an object obtained with the function dist. The function as.dist can be used to transform object of class matrix to object of class dist. -D1 is the dissimilarity matrix D 1 between the same n observations. It must be an object of class dist. By default D1=NULL and the clustering is performed using D0 only. alpha must be a real value between 0 and 1. The mixing parameter α gives the relative importance of D 0 compared to D 1 . By default, this parameter is equal to 0 and only D 0 is used in the clustering process. scale must be a logical value. If TRUE (by default), the dissimilarity matrices D 0 and D 1 are scaled between 0 and 1 (that is divided by their maximum value). wt must be a n-dimensional vector of the weights of the observations. By default, wt=NULL corresponds to the case where all observations are weighted by 1/n.

The function hclustgeo returns an object of class hclust.

The procedure to determine a suitable value for the mixing parameter α is applied through the function choicealpha with the following arguments: choicealpha(D0, D1, range.alpha, K, wt = NULL, scale = TRUE, graph = TRUE)

where:

-D0 is the dissimilarity matrix D 0 of class dist, already defined above.

-D1 is the dissimilarity matrix D 1 of class dist, already defined above.

range.alpha is the vector of the real values α j (between 0 and 1) considered by the user in the grid G of size J.

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

-K is the number of clusters chosen by the user.

wt is the vector of the weights of the n observations, already defined above.

scale is a logical value that allows the user to rescale the dissimilarity matrices D 0 and D 1 , already defined above. graph is a logical value. If graph=TRUE, the two graphics (proportion and normalized proportion of explained inertia) are drawn.

This function returns an object of class choicealpha which contains -Q is a J × 2 real matrix such that the jth row contains Q 0 (P αj K ) and

Q 1 (P αj K ).
-Qnorm is a J × 2 real matrix such that the jth row contains Q * 0 (P αj K ) and Q * 1 (P αj K ).. range.alpha is the vector of the real values α j considered in the G.

A plot method is associated with the class choicealpha.

An illustration of hierarchical clustering with geographical constraints using the package ClustGeo

This section illustrates the procedure of hierarchical clustering with geographical constraints on a real dataset using the package ClustGeo. The complete procedure and methodology for the choice of the mixing parameter α is provided with two types of spatial constraints (with geographical distances and with neighborhood contiguity). We have provided the R code of this case study so that readers can reproduce our methodology and obtain map representations from their own data. employ.rate.city is the employment rate of the municipality, that is the ratio of the number of individuals who have a job to the population of working age (generally defined, for the purposes of international comparison, as persons of between 15 and 64 years of age). graduate.rate refers to the level of education of the population, i.e. the highest qualification declared by the individual. It is defined here as the ratio for the whole population having completed a diploma equal to or greater than two years of higher education (DUT, BTS, DEUG, nursing and social training courses, la licence, maîtrise, DEA, DESS, doctorate, or Grande Ecole diploma). housing.appart is the ratio of apartment housing.

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

agri.land is the part of agricultural area of the municipality.

We consider here two dissimilarity matrices:

-D 0 is the Euclidean distance matrix between the n municipalities performed with the p = 4 available socio-economic variables, -D 1 is a second dissimilarity matrix used to take the geographical proximity between the n municipalities into account.

> # the socio-economic distances > D1 <-as.dist(estuary$D.geo) # the geographic distances between the municipalities

Choice of the number K of clusters

To choose the suitable number K of clusters, we focus on the Ward dendrogram based on the p = 4 socio-economic variables, that is using D 0 only.

> tree <-hclustgeo(D0) > plot(tree,hang=-1,label=FALSE, xlab="",sub="", main="") > rect.hclust(tree,k=5,border=c(4,5,3,2,1)) > legend("topright",legend=paste("cluster",1:5),fill=1:5,bty="n",border="white")

The visual inspection of the dendrogram in Figure 1 suggests to retain K = 5 clusters. We can use the map provided in the estuary data to visualize the corresponding partition in five clusters, called P5 hereafter.

> P5 <-cutree(tree,5) # cut the dendrogram to get the partition in 5 clusters > sp::plot(estuary$map,border="grey",col=P5) # plot an object of class sp > legend("topleft",legend=paste("cluster",1:5),fill=1:5,bty="n",border="white")

Figure 2 shows that municipalities of cluster 5 are geographically compact, corresponding to Bordeaux and the 15 municipalities of its suburban area and Arcachon. On the contrary, municipalities in cluster 3 are scattered over a wider geographical area from North to South of the study area. The composition of each cluster is easily obtained, as shown for cluster 5:

# list of the municipalities in cluster 5 > city_label <-as.vector(estuary$map$"NOM_COMM") > city_label[which(P5==5)] [1] "ARCACHON" "BASSENS" "BEGLES"

Author-produced version of the article published in Computational Statistics, 2018, Fig. 1 Dendrogram of the n = 303 municipalities based on the p = 4 socio-economic variables (that is using D 0 only).

cluster 1 cluster 2 cluster 3 cluster 4 cluster 5

Fig. 2 Map of the partition P5 in K = 5 clusters only based on the socio-economic variables (that is using D 0 only).

The interpretation of the clusters according to the initial socio-economic variables is interesting. Figure 7 shows the boxplots of the variables for each cluster of the partition (left column). Cluster 5 corresponds to urban municipalities, Bordeaux and its outskirts plus Arcachon, with a relatively high graduate rate but low employment rate. Agricultural land is scarce and municipalities have a high proportion of apartments. Cluster 2 corresponds to suburban municipalities (north of Royan; north of Bordeaux close to the Gironde estuary) with mean levels of employment and graduates, a low proportion of apartments, more detached properties, and very high proportions of farmland. Cluster 4 corresponds to municipalities located in the Landes forest. Both the graduate rate and the ratio of the number of individuals in employment are high (greater than the mean value of the study area). The number of apartments is quite low and the agricultural areas are higher to the mean value of the zone. Cluster 1 corresponds to municipalities on the banks of the Gironde estuary. The proportion of farmland is higher than in the other clusters. On the contrary, the number of apartments is the lowest. However this cluster also has both the lowest employment rate and the lowest graduate rate. Cluster 3 is geographically sparse. It has the highest employment rate of the study area, a graduate rate similar to that of cluster 2, and a collective housing rate equivalent to that of cluster 4. The agricultural area is low.

Obtaining a partition taking into account the geographical constraints

To obtain more geographically compact clusters, we can now introduce the matrix D 1 of geographical distances into hclustgeo. This requires a mixing parameter to be selected α to improve the geographical cohesion of the 5 clusters without adversely affecting socio-economic cohesion.

Choice of the mixing parameter α. The mixing parameter α ∈ [0, 1] sets the importance of D 0 and D 1 in the clustering process. When α = 0 the geographical dissimilarities are not taken into account and when α = 1 it is the socio-economic distances which are not taken into account and the clusters are obtained with the geographical distances only. The idea is to perform separate calculations for socio-economic homogeneity and the geographic cohesion of the partitions obtained for a range of different values of α and a given number of clusters K.

To achieve this, we can plot the quality criterion Q 0 and Q 1 of the partitions P α K obtained with different values of α ∈ [0, 1] and choose the value of α which is a trade-off between the lost of socio-economic homogeneity and the gain of geographic cohesion. We use the function choicealpha for this purpose.

> cr <-choicealpha(D0,D1,range.alpha=seq(0,1,0. Figure 3 gives the plot of the proportion of explained pseudo-inertia calculated with D 0 (the socio-economic distances) which is equal to 0.81 when α = 0 and decreases when α increases (black solid line). On the contrary, the proportion of explained pseudo-inertia calculated with D 1 (the geographical distances) is equal to 0.87 when α = 1 and decreases when α decreases (dashed line).

Here, the plot would appear to suggest choosing α = 0.2 which corresponds to a loss of only 7% of socio-economic homogeneity, and a 17% increase in geographical homogeneity.

Final partition obtained with α = 0.2. We perform hclustgeo with D 0 and D 1 and α = 0.2 and cut the tree to get the new partition in five clusters, called P5bis hereafter.

> tree <-hclustgeo(D0,D1,alpha=0.2) > P5bis <-cutree(tree,5) > sp::plot(estuary$map,border="grey",col=P5bis) > legend("topleft",legend=paste("cluster",1:5),fill=1:5,bty="n",border="white")

The increased geographical cohesion of this partition P5bis can be seen in Figure 4. Figure 7 shows the boxplots of the variables for each cluster of the partition P5bis (middle column). Cluster 5 of P5bis is identical to cluster 5 of P5 with the Blaye municipality added in. Cluster 1 keeps the same interpretation as in P5 but has gained spatial homogeneity. It is now clearly located on the banks of the Gironde estuary, especially on the north bank. The same applies for cluster 2 especially for municipalities between Bordeaux and the estuary. Both clusters 3 and 4 have changed significantly. Cluster 3 is now a spatially compact zone, located predominantly in the Médoc. It would appear that these two clusters have been separated based on proportion of farmland, because the municipalities in cluster 3 have above-average proportions of this type of land, while cluster 4 has the lowest proportion of Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1 Two regions with contiguous boundaries, that is sharing one or more boundary point, are considered as neighbors. Let us first build the adjacency matrix A.

> list.nb <-spdep::poly2nb(estuary$map,row.names=rownames(estuary$dat)) #list of neighbors It is possible to obtain the list of the neighbors of a specific city. For instance, the neighbors of Bordeaux (which is the 117th city in the R data table) is given by the script:

> city_label[list.nb [[117]]] # list of the neighbors of BORDEAUX [START_REF] Ambroise | Clustering of spatial data by the EM algorithm[END_REF] "BASSENS" "BEGLES" "BLANQUEFORT" "LE BOUSCAT" "BRUGES" [6] "CENON" "EYSINES" "FLOIRAC" "LORMONT" "MERIGNAC" [START_REF] Legendre | const.clust: Space-and Time-Constrained Clustering Package[END_REF] "PESSAC" "TALENCE"

The dissimilarity matrix D 1 is constructed based on the adjacency matrix A with D 1 = 1 n -A.

> A <-spdep::nb2mat(list.nb,style="B") # build the adjacency matrix

> diag(A) <-1 > colnames(A) <-rownames(A) <-city_label > D1 <-1-A > D1[1:2,1:5] ARCES ARVERT BALANZAC BARZAN BOIS ARCES 0 1 1 0 1 ARVERT 1 0 1 1 1 > D1 <-as.dist(D1)
Choice of the mixing parameter α. The same procedure for the choice of α is then used with this neighborhood dissimilarity matrix D 1 .

> cr <-choicealpha(D0,D1,range.alpha=seq(0,1,0. Fig. 5 Choice of α for a partition in K = 5 clusters when D1 is the neighborhood dissimilarity matrix between municipalities. Top: proportion of explained pseudo-inertias Q 0 (P α K ) versus α (in black solid line) and Q 1 (P α K ) versus α (in dashed line). Bottom: normalized proportion of explained pseudo-inertias Q * 0 (P α K ) versus α (in black solid line) and Q * 1 (P α K ) versus α (in dashed line).

With these kinds of local dissimilarities in D 1 , the neighborhood within-cluster cohesion is always very small. Q 1 (P α K ) takes small values: see the dashed line of Q 1 (P α K ) versus α at the top of Figure 5. To overcome this problem, the user can plot the normalized proportion of explained inertias (that is Q * 0 (P α K ) and Q * 1 (P α K )) instead of the proportion of explained inertias (that is Q 0 (P α K ) and Q 1 (P α K )). At the bottom of Figure 5, the plot of the normalized proportion of explained inertias suggests to retain α = 0.2 or 0.3. The value α = 0.2 slightly favors the socio-economic homogeneity versus the geographical homogeneity. According to the priority given in this application to the socio-economic aspects, the final partition is obtained with α = 0.2.

Final partition obtained with α = 0.2. It remains only to determine this final partition for K = 5 clusters and α = 0.2, called P5ter hereafter. The corresponding map is given in Figure 6.

> tree <-hclustgeo(D0,D1,alpha=0.2)

Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1 > P5ter <-cutree(tree,5) > sp::plot(estuary$map,border="grey",col=P5ter) > legend("topleft",legend=paste("cluster",1:5),fill=1:5,bty="n",border="white") cluster 1 cluster 2 cluster 3 cluster 4 cluster 5

Fig. 6 Map of the partition P5ter in K = 5 clusters based on the socio-economic distances D 0 and the "neighborhood" distances of the municipalities D 1 with α = 0.2.

Figure 6 shows that clusters of P5ter are spatially more compact than that of P5bis. This is not surprising since this approach builds dissimilarities from the adjacency matrix which gives more importance to neighborhoods. However since our approach is based on soft contiguity constraints, municipalities that are not neighbors are allowed to be in the same clusters. This is the case for instance for cluster 4 where some municipalities are located in the north of the estuary whereas most are located in the southern area (corresponding to forest areas). The quality of the partition P5ter is slightly worse than that of partition P5ter according to criterion Q 0 (72.69% versus 75.58%). However the boxplots corresponding to partition P5ter given in Figure 7 (right column) are very similar to those of partition P5bis. These two partitions have thus very close interpretations.

Concluding remarks

In this paper, a Ward-like hierarchical clustering algorithm including soft spacial constraints has been introduced and illustrated on a real dataset. The corresponding approach has been implemented in the R package ClustGeo available on the CRAN. When the observations correspond to geographical units (such as a city or a region), it is then possible to represent the clustering obtained on a map regarding the considered spatial constraints. This Wardlike hierarchical clustering method can also be used in many other contexts where the observations do not correspond to geographical units. In that case, the dissimilarity matrix D 1 associated with the "constraint space?" does not correspond to spatial constraints in its current form.

For instance, the user may have at his/her disposal a first set of data of p 0 variables (e.g. socio-economic items) measured on n individuals on which he/she has made a clustering from the associated dissimilarity (or distance) matrix. This user also has a second data set of p 1 new variables (e.g. environmental items) measured on these same n individuals, on which a dissimilarity matrix D 1 can be calculated. Using the ClusGeo approach, it is possible to take this new information into account to refine the initial clustering without basically disrupting it.

4. 1

 1 The data Data were taken from French population censuses conducted by the National Institute of Statistics and Economic Studies (INSEE). The dataset is an extraction of p = 4 quantitative socio-economic variables for a subsample of n = 303 French municipalities located on the atlantic coast between Royan and Mimizan:

  , 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

Fig. 4

 4 Fig. 4 Map of the partition P5bis in K = 5 clusters based on the socio-economic distances D 0 and the geographical distances between the municipalities D 1 with α = 0.2.

  Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

		1.0					based on D0
							based on D1
		0.8
		0.6
	Q
		0.4
		0.2
		0.0
		0.0	0.2	0.4	0.6	0.8	1.0
				alpha
		of 81%				of 8%
		1.0
		0.8
	Qnorm	0.4 0.6
		0.2
			based on D0
		0.0	based on D1
		0.0	0.2	0.4	0.6	0.8	1.0
				alpha
							1),K=5,graph=TRUE)
	> cr$Q # proportion of explained pseudo-inertia
	> cr$Qnorm # normalized proportion of explained pseudo-inertia

  Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1Fig.7Comparison of the final partitions P5, P5bis and P5ter in terms of variables.Author-produced version of the article published in Computational Statistics, 2018, 33(4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : 10.1007/s00180-018-0791-1

				Partition P5, Cluster 1					Partition P5bis, Cluster 1					Partition P5ter, Cluster 1
	agri.land							agri.land								agri.land
	housing.appart							housing.appart								housing.appart
	graduate.rate							graduate.rate								graduate.rate
	employ.rate.city							employ.rate.city								employ.rate.city
	0	20		40	60	80		0	20		40	60	80		0	20		40	60	80
				Partition P5, Cluster 2					Partition P5bis, Cluster 2					Partition P5ter, Cluster 2
	agri.land							agri.land								agri.land
	housing.appart							housing.appart								housing.appart
	graduate.rate							graduate.rate								graduate.rate
	employ.rate.city							employ.rate.city								employ.rate.city
	0	10	20	30	40	50	60	70	0	20		40		60		80	0	20		40	60	80
				Partition P5, Cluster 3					Partition P5bis, Cluster 3					Partition P5ter, Cluster 3
	agri.land							agri.land								agri.land
	housing.appart							housing.appart								housing.appart
	graduate.rate							graduate.rate								graduate.rate
	employ.rate.city							employ.rate.city								employ.rate.city
	0	20		40	60			0	10	20	30	40	50	60	70	0	10	20	30	40	50	60	70
				Partition P5, Cluster 4					Partition P5bis, Cluster 4					Partition P5ter, Cluster 4
	agri.land							agri.land								agri.land
	housing.appart							housing.appart								housing.appart
	graduate.rate							graduate.rate								graduate.rate
	employ.rate.city							employ.rate.city								employ.rate.city
	0	10		20		30		40	0	20	40		60			0	20	40	60
				Partition P5, Cluster 5					Partition P5bis, Cluster 5					Partition P5ter, Cluster 5
	agri.land							agri.land								agri.land
	housing.appart							housing.appart								housing.appart
	graduate.rate							graduate.rate								graduate.rate
	employ.rate.city							employ.rate.city								employ.rate.city
	0	20		40	60			0	20	40		60			0	20	40	60

http://www.rise-group.org/risem/clusterpy/

http://www.math-evry.cnrs.fr/logiciels/bald

Author-produced version of the article published in Computational Statistics, 2018, 33 (4), 1799-1822 The original publication is available at https://link.springer.com/article/10.1007%2Fs00180-018-0791-1 doi : farmland of the whole partition. Cluster 4 is also different because of the increase in clarity both from a spatial and socio-economic point of view. In addition, it contains the southern half of the study area. The ranges of all variables are also lower in the corresponding boxplots.

Obtaining a partition taking into account the neighborhood constraints

Let us construct a different type of matrix D 1 to take neighbouring municipalities into account when clustering the 303 municipalities.

Author-produced version of the article published in Computational Statistics, 2018, 33 (4)


In a recent paper published in PNAS, Fayet et al. [START_REF] Fayet | Evidence of tool use in a seabird[END_REF] reported scarce field observations of two Atlantic puffins (four years apart) apparently scratching their bodies using sticks, which was interpreted by the authors as evidence of tool use in this species. In a short response, Benjamin Farrar [START_REF] Farrar | Evidence of tool use in a seabird? PsyArXiv, 463hk, ver. 5 recommended and peer-reviewed by Peer Community[END_REF] raises serious concerns about this interpretation and proposes simpler, more parsimonious, mechanisms explaining the observed behaviour: a textbook case of Morgan's canon. In virtually all introductory lectures on animal behaviour, students are advised to exercise caution when interpreting empirical data and weighting alternative explanations. We are sometimes prisoner of our assumptions: our desire of beliefs in advanced cognitive skills in non-human species make us more receptive to facts confirming our preconceptions than to simpler, less exciting, interpretations (a phenomenon known as "confirmation bias" in psychology). We must resist the temptation to accept appealing explanations without enough critical thinking. Our students are thus taught to apply the Lloyd Morgan's canon, a variant of one of the most important heuristics in Science, the principle of parsimony or Occam's razor, rephrased by Morgan [3,page 53] in the context of animal behaviour: "In no case may we interpret an action as the outcome of a higher psychical faculty, if it can be interpreted as the outcome of the exercise of one that stands lower in the psychological scale". In absence of evidence to the contrary, one should postulate the simplest cognitive skill consistent with the observed behaviour. While sometimes criticized from an epistemological point of view [START_REF] Meketa | A critique of the principle of cognitive simplicity in comparative cognition[END_REF][START_REF] Fitzpatrick | Against Morgan's Canon[END_REF][START_REF] Starzak | Interpretations without justification: a general argument against Morgan's Canon[END_REF], it remains an essential and largely accepted framework of animal cognition. It has repeatedly proved to be a useful guide in the minefield of comparative psychology. Classical ethology questions related to the existence of, for instance, meta-cognition [START_REF] Arbilly | Constructive anthropomorphism: a functional evolutionary approach to the study of human-like cognitive mechanisms in animals[END_REF], intentionality or problem solving [START_REF] Taylor | An end to insight? New Caledonian crows can spontaneously solve problems without planning their actions[END_REF] have been convincingly investigated using this principle. Yet, there is a downside to this conservative approach. Blind reference to Morgan's canon may narrow our theoretical thinking about animal cognition [START_REF] Arbilly | Constructive anthropomorphism: a functional evolutionary approach to the study of human-like cognitive mechanisms in animals[END_REF][START_REF] De Waal | Are we smart enough to know how smart animals are?[END_REF]. It could be counter-productive to systematically deny advanced cognitive skills in animals. On the contrary, keeping our mind open to unplanned observations, unexpected discoveries, or serendipity [START_REF] Scheffer | The forgotten half of scientific thinking[END_REF], and being prepared to accept new hypotheses, sometimes fairly remote from the dominant paradigm, may be a fruitful research strategy. To quote Darwin's famous letter to Alfred Wallace: "I am a firm believer, that without speculation there is no good and original observation" [START_REF] Darwin | Letter to A. R. Wallace[END_REF]. Brief notes in specialized scientific journals, or even in grey literature (by enthusiast amateur ornithologists, ichthyologists, or entomologists), constitutes a rich array of anecdotal observations. For instance, Sol et al. [START_REF] Sol | Brain size, innovative propensity and migratory behaviour in temperate Palaearctic birds[END_REF] convincingly compared the innovation propensity across bird species by screening ornithology literature using keywords like 'never reported', 'not seen before', 'first report', 'unusual' or 'novel'. Even if "the plural of anecdote is not data" as the saying goes, such descriptions of novel behaviours, even singlesubject observations, are indisputably precious: taxonomic ubiquity of a behaviour is a powerful argument in favour of evolutionary convergence. Of course, a race to the bottom, amplified by the inevitable media hypes around scientific articles questioning human exceptionalism, is another possible scientific trap for behavioural biologists in search of skills characteristic of so-called advanced species, but never described so far in supposedly cognitively simpler organisms. As stated by Franz de Waal [START_REF] De Waal | Are we smart enough to know how smart animals are?[END_REF]: "I have nothing against anecdotes, especially if they have been caught on camera or come from reputable observers who know their animals; but I do view them as a starting point of research, never an end point". In the case of the two video observations of puffins apparently using sticks as scratching tool, it must be considered as a mere anecdote unless scientists systematically investigate this behaviour. In his constructive criticism of Fayet et al.'s paper, Benjamin Farrar [START_REF] Farrar | Evidence of tool use in a seabird? PsyArXiv, 463hk, ver. 5 recommended and peer-reviewed by Peer Community[END_REF] proposes interesting directions of research and testable predictions. A correlation between the background rate of stick picking and the rate of stick preening would indicate that this behaviour was more likely explained by fluke than genuine innovation in this species.

Benjamin Farrar suggests that the case of tool used reported in puffins and published in PNAS by Fayet and colleagues is not necessarily evidencing tool use in this group of birds. Farrar suggests that an alternative hypothesis to "tool use" could explain this behaviour. Namely, the bird could have picked up a stick for another purpose and coincidentally happen to need to groom himself at the same moment. I fully agree with this more parsimonious hypothesis and with his other point: only 2 observations were documented over several years in this species. I am shocked that a journal like PNAS accepted to publish this work based on so little behavioural elements, and I think that it is really important to publish B Farrar's comments of this study.

Reviewed by Alex Taylor, 2020-01-24 18:08 This is a fair commentary on the interpretation of the tool-preening behaviour recently reported in Fayet et al (2019). I applaud the author for making a testable prediction for their first critique, namely that the background rate of stick pick up could provide useful data for determining whether stick preening was due to chance or not. I would like to see the author do the same for their second critique, that namely testing if the animal had an intention to preen itself with the stick.

Minor comments I think this sentence could do with some unpacking "This should lead us to have very low priors that a puffin would use a stick to scratch itself". Not all readers will know what 'priors' mean in this sense, and so this should be rephrased so it is clearer what point the author is making

  Touchy matter: the delicate balance between Morgan's canon and open-minded description of advanced cognitive skills in the animal Francois-Xavier Dechaume-Moncharmont based on reviews by Valérie Dufour and Alex Taylor A recommendation of: Benjamin G. Farrar. Evidence of tool use in a seabird? (2020), PsyArXiv, 463hk, ver. 5 recommended and peer-reviewed by Peer Community In Ecology. 10.31234/osf.io/463hk Submitted: 22 January 2020, Recommended: 31 January 2020 Cite this recommendation as: Francois-Xavier Dechaume-Moncharmont (2020) Touchy matter: the delicate balance between Morgan's canon and open-minded description of advanced cognitive skills in the animal. Peer Community in Ecology, 100042. 10.24072/pci.ecology.100042

Author's reply:

  Dear Prof. Dechaume-Moncharmont, I would like to thank yourself, Dr Dufour and Dr Taylor for your comments on the article. In-line with Dr Taylor's suggestions I have clarified my use of the word priors, and have made tentative suggestions for how researchers may test whether the puffins were using the sticks with the intention to preen. I have made a couple of suggestions, one involving the continued observation of puffins in the wild and one involving captive animals with access to many sticks, that would open up the possibility for controlled experiments to further test the tooluse hypothesis. Kind regards, Ben Farrar


Lucas Breder Teixeira
email: lucas.breder-teixeira@univ-orleans.fr

Eric Blond

Thomas Sayet

Jean Gillibert

Inverse Identification of Mechanical Materials' Properties Using the Integrated-DIC Technique

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

-Introduction:

• The modelling of refractory materials at high temperatures requires the determination of the most suitable mathematical models

• Moreover, accurate models tend to require complicated testing for the determination of the materials' parameters

• Digital Image Correlation techniques have been used to close the gap between complex materials' models and the determination of the required mechanical properties

-Objectives:

• This work intends to apply the I-DIC technique to the determination of mechanical properties of refractories

• Especial emphasis will be given to their primary and secondary creep behaviour at high temperatures • It's expected to fully characterize the materials (tension, compression and shear) using a reduced number of tests, making use of the Brazilian test

-Conclusions:

• The I-DIC technique is a convenient option to perform inverse identifications, and the robustness of the optimization algorithm related to it plays a decisive role in the obtained results

• For future works, the challenges of the high temperature I-DIC will be addressed, and the development of innovative techniques related to the images' quality is a priority

-I-DIC algorithm:

Convergence?

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement no.764987


Lucas Breder Teixeira
email: lucas.breder-teixeira@univ-orleans.fr

Eric Blond

Thomas Sayet

Jean Gillibert

Inverse Identification of Asymmetric Creep Parameters Using the I-DIC Technique

-Introduction:

• The creep characterization of refractory materials at high temperatures requires complex tests where 𝜎 is the yield stress, 𝒔 is the deviatoric stress and 𝜎 is the von Mises stress. A Genetic Algorithm was used to perform the inverse identification of the parameters. Tables 1 and2 show the reference values used to create the deformed images, the range of variation used in the optimization and the identification results.

-Results:

2 -Brazilian test:

F Shear Tension Compression
The Brazilian test presents a complex stress distribution, with tension, compression and shear observed in the sample at the same time. This allows the identification of the asymmetric behavior using a reduce number of tests

•

  This work intends to apply the Integrated Digital Image Correlation technique to the determination of the secondary creep properties of refractories at high temperatures • It's expected to fully characterize the materials (tension, compression and shear) using a reduced number of tests, making use of the Brazilian test • Images taken at 1200°C were virtually deformed using an asymmetric creep law (different behavior under tension and compression) with known materials' parameters. Later, starting from an initial guess for the parameters, an inverse identification was run to check the I-DIC software robustness 7 -Conclusions: • The I-DIC technique is a convenient option to perform inverse identification of mechanical parameters at high temperatures • For future works, the creep parameters of refractory materials will be identified using experimental images 3 High temperature test setup and Brazilian test samples used with I-DIC. (b) Virtual image. (c) Room temperature. (d) 1200°C without using the blue lights. (e) 1200°C using the blue lights. (f) Images' histograms.

  𝐾 , 𝑛 , 𝑛 )

Table 2 -

 2 Three best solutions given by the Genetic Algorithm

Table 1 -

 1 Reference values for the creep properties Reference valuesRange of variation

E. Blond, N. Schmitt, F. Hild, P. Blumenfeld, J. Poirier, Journal of the European Ceramic Society,

25[11] 1819-1827 (2005).

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement no.764987