Nomenclature

• The term “objects” is used here to refer to the objects which are to be aggregated into clusters; these may be points, lines, polygons, or any other spatial or non-spatial entities.
• The term “non-spatial data” refers to the data on which clusters are based. These data are not necessarily devoid of spatial structure, but are distinct from the spatial coordinates used to constrain the clustering.

Spatial clustering

Spatial clustering is a very well-studied field (see reviews in Han et al. 2001; Duque et al. 2007; Lu 2009), with many methods implemented in the R language (see the CRAN Task View on Analysis of Spatial Data). Spatial clustering algorithms take as input a set of spatial distances between objects, and seek to cluster those objects based on these exclusively spatial distances alone (Fig. 1A). Other non-spatial data may be included, but must be somehow reconciled with the spatial component. This is often achieved through weighted addition to attain approximately “spatialized” data. Figure 1B-C illustrate two related non-spatial (B) and spatial (C) datasets. The primary data of interest depicted in B can be “spatialized” through additively combining the associated distance matrices of non-spatial and spatial data, and submitting the resultant distance matrix to a clustering routine of choice.

Figure 1: (A) Illustration of typical spatial clustering application for which input data are explicit spatial distances between points; (B) Illustration of clustering in some other, non-spatial dimensions, D1 and D2, for which associated spatial data in (C) do not manifest clear spatial clusters.

For example, the following code illustrates the use of the DBSCAN algorithm (Density Based Spatial Clustering of Applications with Noise) from the R package dbscan.

d_nospace <- dist (dat_nospace) # matrix of non-spatial data
d_space <- dist (dat_space) # 2-column matrix of spatial data
d <- d_nospace + d_space # simple linear addition

library (dbscan)
db <- dbscan::dbscan (d, eps = 0.4) # more on the eps parameter below
db

## DBSCAN clustering for 49 objects.
## Parameters: eps = 0.4, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 4 cluster(s) and 4 noise points.
## 
##  0  1  2  3  4 
##  4 17  9  7 12 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints

Figure 2: (A) Non-spatial data coloured by dbscan clustering results; (B) Corresponding spatial data coloured by dbscan clustering results.

These results demonstrate that reasonable results can indeed be obtained through simple linear combination of non-spatial and spatial distances. This approach is very simple, and it is very easy to submit such combined distance matrices to high-performance clustering routines such as dbscan. There are nevertheless two notable shortcomings:

1. There is no objectively best way to combine non-spatial and spatial distance matrices; and
2. Routines such as dbscan still require an effectively arbitrary parameter represented by the above value of eps = 0.4. This value was simply chosen to best reflect the known structure of the input clusters, but in any practical application will remain largely arbitrary.

Even one of the most recent R packages dedicated to spatial clustering (Chavent et al. 2017) follows precisely this linear addition strategy, via a parameter determining the relative weights of the non-spatial and spatial data.

We assert here that such approaches provide more a means of attaining approximately spatially-structured clustering schemes, rather than actually providing spatially-constrained clusters in the sense we now explore.

Spatially-constrained clustering

The spatialcluster package performs strict spatially constrained clustering. This means that clusters are formed entirely on the basis of the non-spatial data, while the spatial data provide a constraint used to ensure that all clusters are spatially contiguous.

The SKATER algorithm

One of the most widespread algorithms for spatially constrained clustering is the so-called “SKATER” algorithm (Spatial Kluster Analysis by Tree Edge Removal, AssunÇão et al. 2006), available via the R package spdep (Bivand et al. 2013; Bivand and Piras 2015). This algorithm constructs a minimum spanning tree (MST) connecting neighbours defined by the spatial data, with all connections defined by minimal distances in the non-spatial data. The SKATER algorithm constructs the simplest of all MSTs, by iterating through the list of neighbouring edges according to increasing (non-spatial) distance, and inserting edges where these connect previously unconnected objects. The resultant MST is then partitioned into a specified number of clusters such that the intra-cluster sum of squared deviations from the mean (of the non-spatial data) is minimised. Further details are given in AssunÇão et al. (2006).

The REDCAP algorithms

The REDCAP algorithms for spatially-constrained clustering (REgionalization with Dynamically Constrained Agglomerative Clustering and Partitioning, Guo 2008) employ three distinct methods for constructing MSTs. The original algorithms develop these three methods for two distinct ways of constructing spanning trees: through using nearest neighbours only (“first-order constraint”), or through considering all neighbours of each object (“full-order constraint”). The results of the cited manuscript clearly reveal the superiority of the latter, which is the default here, although first-order constraints may be requested via the full_order parameter of scl_redcap().

The three methods refer to methods for determining which edges are selected to link clusters, via either single-, average-, or complete-linkage, which function as follows:

1. Single-linkage clustering simply loops through the list of minimal-distance, nearest-neighbour edges, and inserts each next-minima distance edge into the clustering tree if it is not part of any previous cluster, and if it connects two separate yet contiguous clusters.
2. Average-linkage clustering assigns weights to each unassigned edge based on their average distance to all edges within all adjacent clusters. Thus when an edge becomes part of a cluster, the distances to all non-assigned edges adjacent to that cluster are updated to reflect the change in average distance calculated over all edges in that cluster. Edges are continually (re-)sorted based on average distances, and the tree is built through sequentially inserting minimal-distance edges.
3. Complete-linkage (or maximal-linkage) clustering forms clusters through inserting the edge having the minimal distance to the farthest point of any adjacent cluster.

Single-linkage is equivalent to the SKATER algorithm, where the single best edge connecting two clusters is selected. (The SKATER algorithm is actually equivalent to the worst-performing REDCAP algorithm: single-linkage, first-order constraint.) See Guo (2008) for details. The resultant MSTs are then partitioned into specified numbers of clusters using an identical approach to SKATER, namely through minimising the intra-cluster sum of squared deviations from mean values.

The “full” clustering algorithm

As described above, the REDCAP algorithms provide different ways of constructing minimal spanning trees for a given data set, with resultant clusters based on bisecting these spanning trees. They effectively represent the relationships within a data set by the “best” minimal set (according to the chosen algorithm). The package also implements a simple alternative, here called the “full” algorithm, for deriving a clustering scheme that uses the full set of nearest-neighbour relationships in a given data set.

Each point or object within a typical planar (effectively two-dimensional) data set may have on average just under three nearest neighbours if these are calculated with a triangulation, or potentially up to \(k\) for some \(k\)-nearest neighbours scheme. Each point in a minimal-spanning tree generally has an average of between two and three neighbours (one for terminal nodes; two for non-branching nodes; three for branching nodes). It may accordingly be presumed that reducing a full set of neighbours to an MST reduces average numbers of edges per node from \(\gtrapprox 3\) to \(\sim 2.5\). The loss in computational complexity produced through using all neighbouring nodes instead of an MST is thus likely to be only around 20%. Moreover, if a clustering algorithm scales sub-linearly with \(N\), as many do, this loss is likely to be even less pronounced with increasing \(N\).

The algorithm proceeds through the following steps, looping until all data objects have been allocated to a cluster, and beginning with both i = 0 and clnum = 0. The primary data is an array of edges (edge) sorted by increasing distance.

1.  Select edge [i] connecting nodes a and b.
2.  i++
3.  if both a and b are already in clusters:
        continue to next iteration
    else if neither node is in a cluster:
        cluster (clnum) = create_new_cluster (c (a, b))
        clnum++
    else if only a is in a cluster:
        set_cluster_number (a, get_cluster_number (b))
    else if only b is in a cluster:
        set_cluster_number (b, get_cluster_number (a))

This procedure allocates all nodes (data objects) to clusters. The number of resultant clusters can not be known in advance, and many clusters may be very small, comprising for example only a single edge connecting two nodes. This initial clustering is then used as the basis of a further hierarchical clustering, through sequentially merging clusters having the minimal distance according to some specified metric. This merging involves linking previously distinct clusters according to some specified criteria, for which the same three criteria used in the REDCAP algorithms can be directly applied here. Note, however, that average and maximal distances in the REDCAP algorithms simply represent the respective average and maximal individual edge distances, and not the corresponding distances traversed within a cluster through the MST.

It would be possible to select edges based on actual average or maximal traversal distances from a potential connecting edge to all other edges in two clusters. Doing so, however, simply merges the two smallest clusters, all other things being equal, because these must by definition have the shortest average and maximal distance. Edge selection is therefore implemented in the same way as the REDCAP algorithms, by selecting edges based on average and maximal single edge distances. Clusters are hierarchically connected by selecting edges according to one of the following three schemes:

1. Single-linkage: select the single edge having the minimal distance between any two clusters;
2. Average-linkage: select the edge connecting the two clusters which, when merged, give the lowest average intra-cluster edge distance; or
3. Complete-linkage: select the edge connecting the two clusters which, when merged, gives the lowest maximal intra-cluster edge distance.

The latter two of these yield clusters with preferentially shorter intra-cluster distances, yet the selection procedure remains statistically unbiased by cluster size.

Using the package

The package provides two main clustering functions, scl_redcap() for the REDCAP algorithms (and thus SKATER), and scl_full() for the full algorithm described above. Both take a two-column matrix of spatial coordinates, xy, a square matrix, dmat, of the non-spatial data to be clustered, and a desired number of clusters, ncl. They return an object of class scl, which can be plotted with the generic plot() method. The following code illustrates usage with random data.

set.seed (1)
n <- 100
xy <- matrix (runif (2 * n), ncol = 2)
dmat <- matrix (runif (n^2), ncol = n)
scl <- scl_full (xy, dmat, ncl = 6)
plot (scl)

Because both algorithms build a complete hierarchy of merges, the same object can be cut into a smaller number of clusters without repeating the computationally expensive construction of the tree, using scl_recluster(); and the hierarchy of merges itself can be visualised with plot_merges().

scl <- scl_recluster (scl, ncl = 3)
plot (scl)

Observation and reference matrices for covariance clustering

The motivating application is spatially-constrained clustering based on non-spatial data representing objects which also have spatial locations. “Spatially-constrained” implies that the clusters are constrained to be spatially contiguous. Spatial data that involve counts or densities are often represented in terms of flow matrices or, in conventional studies of human transport systems, “origin-destination” matrices. Such matrices are rectangular, and represent tallies of numbers, densities, or, in general terms, flows between a set of origins and a (possibly different) set of destinations.

Such flow data can - and indeed often are - directly submitted to clustering routines, through constructing some suitably-scaled metric of distance as the inverse of flow. In spatial systems, however, flows from some point \(A\) to another point \(B\) are often closely related to flows from some point \(A^\prime\) near \(A\) to some \(B^\prime\) near \(B\), with similarity positively related to the proximity of \(A^\prime\) to \(A\) and \(B^\prime\) to \(B\). This is nothing other than a definition of “autocorrelation”. Clustering with directly observed flows alone ignores the effects of autocorrelation, and thus introduces a host of well-recognised issues arising through doing so.

Flow matrices may be converted into covariance matrices through populating the upper or lower triangle with pair-wise covariances (or through copying one triangle to the other to generate a symmetrical covariance matrix). Covariances can be calculated between the rows of a flow matrix to yield covariances in flows from each point, or between the columns to yield covariances in flows to each point. Both may be combined through, for example, filling the upper triangle of a covariance matrix with row-wise covariances, and the lower triangle with column-wise covariances. The value \(C_{ij}\) will then represent the covariance in flow from i and flow from j, while \(C_{ji}\) will represent the covariance in flow to i and flow to j. Such compound covariance matrices will not be symmetrical.

Values of covariance will manifest similar properties to raw flows, in that nearby values are likely to be more similar or “autocorrelated.” For many spatial processes, “neutral” models of flows have been developed, and can be applied to yield estimates of covariances expected in the presence of these kinds of processes alone. Such estimated covariance matrices are the “reference” matrices referred to above, and are generated here through applying spatial interaction models to the observed flow data, and generating corresponding covariance matrices. These matrices effectively represent the covariance expected in the presence of a particular type of autocorrelation. The difference with observed covariance matrices may then be interpreted to reflect the effects of processes beyond spatial autocorrelation.

Observation and reference matrices for correlation clustering

The ability to calculate observation and reference covariance matrices enables a standardised covariance matrix to be generated, in which positive values represent values of covariance beyond those expected on the basis of the neutral model alone (here, a spatial interaction model). Note, moreover, that negative values may also be interpreted in the obverse way: to reflect spatial clusters in which flows are significantly less likely than those expected on the basis of a neutral model alone. The methods described here are applied separately to positive and negative values of the net covariance matrix, with clusters only formed from either exclusively positive or exclusively negative values. Covariances in such clusters will thus always be entirely significant in the sense of T-tests for divergence from expected values of zero. Thus by definition, any clusters discerned by the methodology may be presumed entirely significant. The software produces summary statistics for each cluster (T- and p-values), enabling the relative significance of different clusters to be compared.

Note that it will not generally be possible - absent specific hypotheses - to provide a global statistic quantifying the significance of a clustering scheme. This is primarily because there is no reason to presume that positive net covariances will be more likely to form distinct clusters than will negative net covariances. The overall distribution of net covariances will thus generally be symmetrical, although it may also be asymmetrical when spatial patterns of positive net covariance differ from those of negative net covariance - for example through positive covariances being very strongly concentrated in a small number of clusters, while negative covariances remain more dispersed.

Interpretation of OD clusters

As described above, clusters discerned from net covariance matrices derived from flow data may be either significantly positive or significantly negative. The results of such analyses must be interpreted with caution. Although all clusters may by definition be considered individually significant, these are clusters of covariances and not of underlying flows. Thus, for example, positive clusters calculated from origin data (row-wise covariances) will not reflect regions from which movement is significantly more likely. Rather, each point in such a cluster will be a point from which movement covaries with all other clustered points more than expected on the basis of a neutral model. A positive cluster calculated from origin data may thus be interpreted to represent a distinct region from which movement significantly covaries in a manner independent of the precise point of origin. Importantly, the significance of these net origin- or destination-clusters does not reflect any net patterns of in- or out-flow: net origins (destinations) may still be regions of net inflow (outflow). Rather, these will be regions towards or away from which patterns of flow or movement covary more closely than neutrally expected, independent of actual net rates of flow or movement.

For clarity, the four kinds of cluster able to be discerned with these analyses are as follows:

1. Positive clusters calculated from origin data represent regions from which movement covaries significantly more than expected, and thus represent significant net origins - and importantly, this phenomenon is independent of net rates of out- or in-flow;
2. Negative clusters calculated from origin data represent regions from which movement covaries significantly less than expected, and thus represent significant net non-origins from which movement is expected to covary yet doesn’t;
3. Positive clusters calculated from destination data represent regions toward which movement covaries significantly more than expected, and thus represent significant net destinations; and
4. Negative clusters calculated from destination data represent regions toward which movement covaries significantly less than expected, and thus represent significant net non-destinations.

The phrases “non-origins” and “non-destinations” refer to regions within which movement respectively from or to adjacent locations covaries significantly less than expected. Hypothetical causes are likely to be highly specific for any given study, but may include for example the dominance of a single location within these clusters to or from which movement is focussed at the expense of neighbouring locations. Such dominant locations would then represent locations which effectively act to weaken general patterns of spatial autocorrelation (or interaction) which would otherwise exist.

A formal argument

The preceding intuition can be made precise by considering destinations only; the argument for origins is identical, with rows replacing columns. Let \(F\) be an \(n \times m\) observed flow matrix with columns \(f_a \in \mathbb{R}^n\), where \(f_{ia}\) is the flow from origin \(i\) to destination \(a\), and let \(S\) be the corresponding spatial interaction model with columns \(s_a\). Write \(\mathbf{1}\) for the vector of \(n\) ones, and \(\bar{f}_a = \mathbf{1}^\top f_a / n\) for the mean flow to \(a\).

A net-flow measure sees only one dimension. Any measure of net flow beyond that expected, whether absolute or relative, is a function of the aggregate \(\delta_a = \bar{f}_a - \bar{s}_a = \mathbf{1}^\top (f_a - s_a) / n\). This is the projection of the residual vector \(f_a - s_a \in \mathbb{R}^n\) onto the single direction \(\mathbf{1}\). Two destinations \(a\) and \(a^\prime\) are therefore indistinguishable to such a measure whenever their residuals have equal sums, however different the vectors themselves may be.

A covariance sees the remaining \(n - 1\) dimensions. The covariance between destinations is \[ C_{aa^\prime} = \frac{1}{n} \sum_i (f_{ia} - \bar{f}_a)(f_{ia^\prime} - \bar{f}_{a^\prime}) = \frac{1}{n} \langle P f_a, P f_{a^\prime} \rangle, \qquad P = I - \tfrac{1}{n} \mathbf{1} \mathbf{1}^\top, \] where \(P\) projects onto the orthogonal complement of \(\mathbf{1}\). It thus uses exactly the information discarded by the net-flow measure, and is unchanged by adding any constant to all flows to a destination (so is independent of the net rate of in- or out-flow, as stated above). By the Cauchy-Schwarz inequality, \(|C_{aa^\prime}| \le \sigma_a \sigma_{a^\prime}\), with equality only if \(P f_a\) and \(P f_{a^\prime}\) are parallel; that is, only if flows to \(a\) and \(a^\prime\) are related by an affine function at every origin, and not merely in aggregate.

Spatial contiguity does not remedy this for net flows. Suppose the origins are partitioned into two sub-populations, \(I_1\) and \(I_2\) with \(|I_1| = |I_2| = n/2\), and that destination \(a\) attracts flow \(x\) only from \(I_1\), while neighbouring \(a^\prime\) attracts the same flow \(x\) only from \(I_2\). A spatial interaction model depends only on the sizes of, and distances between, origins and destinations, and assigns the same flows \(s_a = s_{a^\prime}\) to the two neighbours. Then \[ \bar{f}_a = \bar{f}_{a^\prime} = x / 2, \quad \text{so} \quad \delta_a = \delta_{a^\prime}, \] and the two destinations are identical in terms of net flow, and would be merged first by any of the linkage schemes above. Yet the supports of \(f_a\) and \(f_{a^\prime}\) are disjoint, so that \(\sum_i f_{ia} f_{ia^\prime} = 0\) and \[ C_{aa^\prime} = 0 - \bar{f}_a \bar{f}_{a^\prime} = -\frac{x^2}{4} < 0. \] Their flows are in fact maximally anti-correlated, so the two destinations could never be clustered together on the basis of positive net covariance. The spatial constraint determines only which pairs may be merged; it does not ensure that the merged objects share any structure beyond that measured by the chosen dissimilarity. For net flows, that structure is a single scalar, and spatially adjacent objects with equal totals but unrelated origins will be clustered. For covariances, a cluster of mutually positive values \(C_{aa^\prime} > 0\) implies that, on average, the entire population varies in the same direction relative to its mean in its flows to every member of the cluster, which is precisely the interpretation of covariance clusters developed above.