Objectives

• Introduce contiguity constraints into k-means clusering

• Identify contiguous clusters by means of the skater algorithm

• Identify contiguous clusters with the number of clusters as endogenous with max-p

• Comparing the characteristics of the clusters

Getting started

With GeoDa launched and all previous projects closed, we again load the Guerry sample data set from the Connect to Data Source interface. We either load it from the sample data collection and then save the file in a working directory, or we use a previously saved version. The process should yield the familiar themeless base map, showing the 85 French departments.

To carry out the spatially constrained cluster analysis, we will need a spatial weights file, either created from scratch, or loaded from a previous analysis (ideally, contained in a project file). The Weights Manager should have at least one spatial weights file included, e.g., Guerry_85_q for first order queen contiguity, as shown here.

Principle

As we have seen in a previous chapter, the K Means cluster procedure only deals with attribute similarity and does not guarantee that the resulting clusters are spatially contiguous. A number of ad hoc procedures have been suggested to create clusters where all members of the cluster are contiguous, but these are generally unsatisfactory.

An alternative is to include the geometric centroids of the observations as part of the optimization process. There are two general approaches to carry this out. In one, the x and y coordinates are simply added as additional variables in the collection of attributes. While this tends to yield more geographically compact clusters, it does not guarantee contiguity. The second approach consists of a weighted multi-objective optimization, that combines the objective of attribute similarity with that of geometric similarity. This allows for the identification of a cluster with maximal between-cluster dissimilarity that consists of spatially contiguous units.

We consider each approach in turn.

Geometric centroids as attributes

Weighted optimization

The previous approach included the x and y coordinates as individual variables treated in exactly the same manner as the other attributes (the algorithm does not know that these are coordinates). In contrast, in the weighted optimization, the coordinate variables are treated separately from the regular attributes, in the sense that there are now two objective functions. One is focused on the similarity of the regular attributes (in our example, the six variables), the other on the similarity of the geometric centroids. A weight changes the relative importance of each objective.

To illustrate this, we first consider the case where all the weight is given to the geometric coordinates. We select the same six variables as before, not including the coordinates. In the K-Means Settings dialog, we check the box next to Use geometric centroids. For now, we leave the weight to the value of 1.00. This means that the actual attributes (the six variables) are excluded from the optimization process. The K-Means algorithm is applied to the x, y coordinates only. By setting the weight, the coordinates do not need to be specified explicitly, all the calculations are done in background.

The resulting cluster map shows four compact regions with all members of each cluster contiguous to another member.

Even though the clustering algorithm ignores the actual attributes, the resulting properties can still be calculated. As is to be expected, the results are not that good relative to the pure K-Means. The overall ratio of between to total sum of squares is 0.246, compared to 0.433 in the base case.

With the weight set to 0, the same result as in the base case is obtained (not illustrated here). Before moving to an optimization strategy, we check the results with the weight set at 0.5.

The cluster map shows four contiguous subregions, but different from the result obtained in the pure centroid solution.

Relative to the pure centroid solution, the overall characteristics have improved substantially, to 0.347 for the between to total sum of squares ratio, but still well below the benchmark a-spatial k-means.

An optimization bisection search

The selection of a value for the weight in the geometric centroids setting constitutes a compromise between attribute similarity (best solution for w=0) and locational similarity (best solution for w=1). We could thus approach the selection of a value for the weight as an optimization problem, where we try to maximize the between to total ratio of the sum of squares, while ensuring contiguity among cluster members.

In a small example like the 85 French departments, we can pursue this manually by means of a bisection search. The logic behind the bisection search is to enclose the optimal solution between an upper and a lower bound for a parameter. At each step, a mid-point is selected that brings us closer to an optimal solution.

In the cluster example, the natural bounds are 0 (pure K-means) and 1 (pure centroids). The midpoint is 0.50. As we just observed, with this value for the weight, the solution is contiguous. This implies that we can do better on the ratio objective by moving closer to the full K-means solution. We thus take as the next value for the weight the midpoint between 0 and 0.50, i.e., 0.25, and check for contiguity.

The result has a better ratio of 0.417, but two of the regions no longer have contiguous members, with group 4 consisting of four different sub-regions.

Given that 0.25 does not satisfy the contiguity criterion, we now have to move to the right, to the midpoint between 0.25 and 0.50, or 0.375.

This yields a map that almost meets the contiguity constraint, with a ratio of 0.368. Cluster 1 has one observation that is geographically separated from the rest of the cluster.

This means that the next step should be the midpoint between 0.375 and 0.5, or 0.4375. We continue the process until a contiguous solution provides no further meaningful improvement in the sum of squares ratio. In our example, this is for a weight of 0.4500, which yields the following map, with a sum or squares ratio of 0.361.

This is the best we can do in terms of the between to total sum of squares ratio while ensuring contiguity among the members of each cluster. We can now repeat this process for different values of k to assess which combination of number of clusters and centroid weights yields the most satisfactory solution.

Principle

An approach that explicitly takes into account the contiguity constraints in the clustering process is based on the so-called skater algorithm outlined in Assuncao et al. (2006). The algorithm carries out a pruning of the minimum spanning tree created from the spatial weights matrix for the observations. The weights in the spatial weights matrix correspond to the pair-wise dissimilarity between observations, which become the edges in the graph representation of the weights (the observations are the nodes).

Starting from this weights matrix, a minimum spanning tree is obtained, which is a path that connects all observations (nodes), but each is only visited once. In other words, the n nodes are connected by n-1 edges, such that the between-node dissimilarity is minimized. This tree is then pruned by selecting the edge whose removal increases the objective function (between group dissimilarity) the most. The process continues until convergence. This is a hierarchical process, in that once the tree is cut at one point, all subsequent cuts are limited to the resulting subtrees. In other words, once an observation ends up in a pruned branch of the tree, it cannot switch back to a previous branch. This is sometimes viewed as a limitation of this algorithm.

Implementation

In GeoDa, the skater algorithm is invoked from the Clusters toolbar icon, or from the main menu as Clusters > Skater.

Selecting this option brings up the Skater Settings dialog, which has essentially the same structure as the corresponding dialog for K-Means and other classic clustering approaches. The left panel provides a way to select the variables, the number of clusters and different options to determine the clusters. Different from previous approaches is the presence of a Weigths drop down list, where the contiguity weights must be specified. The skater algorithm does not work without a spatial weights file.

The Distance Function, Transformation and seed otpions work in the same manner as for classic clustering. Again, at the bottom of the dialog, we specify the Field or variable name where the classification will be saved.

We consider the Minimum Bound and Min Region Size options below.

We select the same six variables as before: Crm_prs, Crm_prp, Litercy, Donatns, Infants, and Suicids. The spatial weights are based on queen contiguity (Guerry_85_q).

Clicking on Run brings up the skater cluster map as a separate window and lists the cluster characteristics in the Summary panel.

The spatial clusters generated by skater follow the minimum spanning tree closely and tend to reflect the hierarchical nature of the pruning. In this example, the country initially gets divided into three largely horizontal regions, of which the top one gets split into an eastern and a western part.

The ratio of between sum of squares to total sum of squares is 0.316, somewhat lower than the 0.361 we obtained in the bisection weighted optimization above.

The hierarchical nature of the skater algorithm can be further illustrated by increasing the number of clusters to 6 (all other options remain the same). The resulting cluster map shows how the previous regions 2 and 3 each get split into an eastern and a western part.

The ratio of between sum or squares to total sum of squares is now 0.420.

Principle

The max-p region problem, outlined in Duque, Anselin, and Rey (2012), makes the number of regions (k) part of the solution process. This is accomplished by introducing a minimum size constraint for each cluster. In contrast to the skater method, where this was optional, for max-p the constraint is required.

The algorithm itself consists of a search process that starts with an initial feasible solution and iteratively improves upon it while maintaining contiguity among the elements of each cluster. This can take some time, especially for larger data sets. In addition, the iterative process tends to be quite sensitive to the starting point, so an extensive sensitivity analysis is highly recommended.

Implementation

The max-p option is available from the main menu as Clusters > Max-p or as a similarly named option from the Clusters toolbar icon.

The variable settings dialog is similar to the one for the skater algorithm. We select the same six variables. Also, spatial weights need to be specified in the Weights drop down list (here, we take first order queen contiguity as in Guerry_q_85). In addition, we have to set the MIninum Bound variable and threshold. In our example, we again select Pop1831 and take the 10% default. Finally, we set the number of Iterations to 1000 and leave all other options to their default settings. We specify the Field for the resulting cluster membership as CL_p1k, as shown below.

The resulting cluster map yields eight clusters, consisting of 14 to 8 elements.

The corresponding ratio of between sum of squares to total sum of squares is 0.511. The within sum of squares are fairly evenly distributed, ranging from a low of 11.2 for cluster 5 (9 elements) to a high of 49.2 for cluster 6 (also with 9 elements).

Sensitivity analysis

The max-p optimization algorithm is an iterative process, that moves from an initial feasible solution to a superior solution. However, this process may be slow and can get trapped in local optima. Hence, it behooves us to carry out an extensive sensitivity analysis. There are two main approaches to implement this. One is to increase the number of iterations, starting from the same initial feasible solution. An alternative approach is to adjust the initial feasible solution (the starting point for the search), for example, to start from a solution from an earlier set of iterations. Subsequent iterations may improve on this solution.

The motivation behind the sensitivity analysis is to continue iterating until the resulting solution is stable, i.e., it no longer changes and it reaches a maximum for the objective function.