In this document we give a high-level tutorial for using the functionality available in the cordillera package. The technical details are described in Rusch, Hornik & Mair (2018).
The idea of the OPTICS cordillera (OC) is to quantify the “clusteredness” that is inherent in a data matrix \(X\) in \(R^N\). We understand clusteredness to be a spatial property of the overall arrangement of vectors in the \(N\)-dimensional space. There exists a continuum ranging from no clusteredness to perfect clusteredness and we qunatify where on that continuum a data matrix falls. Crucially, it is not the same thing as actually having or assigning clusters, it is more of a structural property of the whole data matrix. To make that distinction clearer, we also refer to this as the “accumulation tendency” and the clusters of points as accumulation. In that it shares some similarity to measures of spatial arrangement, like Ripley’s k-function, only that we try to specifically capture clusteredness. A perhaps more familiar analogy to clusteredness would be whether the vectors in a data matrix are somehow associated or not; the specific type of association being equivalent to how clustering relates to clusteredness. Clusteredness as measured by the OC can be thought of as a summary of all the clustering information in the data matrix and in that being akin to what hierarchical clustering does, only we aggregate this into a single unidimensional measure. If a data matrix has high clusteredness then the clustering can be visually appreciated in a plot and a clustering method applied to that data matrix will typically work well.
The OC is a nonparametric measure in the sense of trying to make as little assumptions as possible. The underlying “cluster” definition is density-based and derived from Ester et al. (1996). It is also underlying the DBSCAN algorithm (Ester et al., 1996). Essentially clusters are defined as accumulations of points that are close to each other (“density-connected” on some distance metric) or as regions with a high density of points that are separated from each other by regions of low or no point density. This definition also allows for nested clusters as regions of higher density within regions of lower density. Furthermore it contains a concept of “noise” points which are points in regions where the overall density is too low for this to count as a cluster. To achieve this only two assumptions must be made: First, that a cluster comprises at least \(k\) observations (minpts) and that the density is assessed within a maximum radius \(\epsilon\) around each point (anything not matching that will count as noise). Due to this setup there is no need to make decisions a priori about the number of accumulations or clusters, nor about the specific shape, intracluster variance or centroid. This illustrates why the concept is appealing for measuring the “clusteredness” as opposed to a specific concrete clustering.
OPTICS
The OPTICS algorithm (Ankerst et al., 1999) is at the heart of the OPTICS cordillera. An R native OPTICS implementation is available in the dbscan package (Hahsler et al., 2019). We also provide a rudimentary interface to the OPTICS reference impementation in ELKI (which would have to be installed for it to work). OPTICS stands for Ordering Points To Infer Clustering Structure and is essentially a hierarchical clustering version of DBSCAN. It abstracts from a concrete \(\varepsilon\) by allowing for a maximum \(\varepsilon\), \(\epsilon\) up to which all \(\varepsilon\) are taken into account (so it just needs to be “large enough”). If OPTICS would be limited to only a specific \(\varepsilon\) (like cutting a dendrogram at a specific height) it would result in the DBSCAN clustering.
Due to this being a hierarchical clustering idea, OPTICS doesn’t give a concrete clustering but orders the points based on a quantity derived from \(k\) and \(\epsilon\) called “reachability”. This ordering is quite complex and happens algorithmically, so it can’t be expressed in close form. What results is the tuple of ordering of the points/row vectors \(x_i\) in \(X\), \(R(X)\), together with the associated reachability \(r_i\) for \(x_i\). This tuple contains the complete clustering structure within \(X\) up to the maximum radius \(\epsilon\) given an accumulation must comparise at least \(k\) points. It holds that points that are subsequent in the ordering and have relatively low reachability belong to the same accumulation, whereas points that are far away from each other in the ordering or subsequent points with relatively high reachabilility belong to different accumulations. Every high reachability signals the beginning of a new accumulation.
This can be plotted with the ordering of the \(x_i\) on the \(x\)-axis and their associated reachability \(r_i\) on the \(y\)-axis. This plot is coined the “reachability plot” and is the visual expression of the complete clustering structure within \(X\) up to the maximum radius \(\epsilon\) given an accumulation must comparise at least \(k\) points. The reachability plot is similar to a dendrogram in that it can be cut at some level \(\varepsilon \leq \epsilon\) to yield a concrete clustering (the corresponding DBSCAN clustering for \(k\) and \(\varepsilon\)).
To illustrate let’s first create a strange cluster situation in 2D.
library(clusterSim)
#> Loading required package: cluster
#> Loading required package: MASS
#three spherical clusters with different variance with n=66 each
n<-198
set.seed(1)
x <- cbind(
x = c(runif(1, -2, 5) + rnorm(n/3, sd=0.2), runif(1, -2, 5) + rnorm(n/3, sd=0.3), runif(1, -2, 5) + rnorm(n/3, sd=0.05)),
y = c(runif(1, -3, 5) + rnorm(n/3, sd=0.2), runif(1, -3, 5) + rnorm(n/3, sd=0.3), runif(1, 0, 4) + rnorm(n/3, sd=0.05))
)
cl<-c(rep(1:3,each=n/3)) #the clusters
#two worm clusters of size 100 each
set.seed(1)
sw<-shapes.worms(100,shape1x1=-2,shape1x2=0.5,shape2x1=-0.3,shape2x2=3)
#noise
set.seed(5)
ns<-cbind(runif(10,-2,5),runif(10,-2,5))
x<-rbind(x,sw$data,ns)
#x<-rbind(x,sw$data)
cl<-c(cl,sw$cluster+3,rep(8,10))
#cl<-c(cl,sw$cluster+3)
plot(x, col=cl,pch=19,xlim=c(-3,5),ylim=c(-3,6))
We see that there are 5 clusters, two are the worms, three are spherical with different variance. The green spherical cluster is nested within the blue worm. Additionally, we have 10 noise points (in grey).
The clustering structure is clear when looking at the plot: The black, red, turquoise, green and blue clusters have a higher density of points than the surrounding area. All clusters but the green cluster are separated by regions of low or no point density. Low density because there are some noise points. The green cluster is a high density region within the blue cluster. This situation is extremely difficult to disentangle for a cluster algorithm that uses centroids, assumes constant variances or spherical shapes. Yet, we see that is quite clustered as an overall plot.
OPTICS can help us with characterizing the situation. Let’s run OPTICS on this, with the minimum number of points comprising a cluster to be \(k=22\) and a large \(\epsilon=10\).
ores<-optics(x,minPts=22,eps=10)
print(ores)
#> OPTICS ordering/clustering for 408 objects.
#> Parameters: minPts = 22, eps = 10, eps_cl = NA, xi = NA
#> Available fields: order, reachdist, coredist, predecessor, minPts, eps,
#> eps_cl, xi
The OPTICS ordering is available as $order and the reachabilities via $reachdist. We can plot those and color the plot based on the known clustering.
The interpretation is now this: Each “valley” in the plot stands for a cluster and each cluster is separated by a peak. Valleys after low peaks are typically nested clusters. The deeper the valley, denser the cluster and the higher the peak, the more separated the clusters. From left to right we see that, the first valley is the black cluster, and then we have a peak signalling the start of the next cluster (the blue one). Then there is a peak signalling a new cluster (the green cluster) but that peak is relatively small compared to the other peaks, so we can assume that the green cluster is nested within the blue cluster (although that is certainly not conclusive here, it could also be that the green cluster is just close to the blue cluster). The next peak is signalling the turqoise cluster and then we have the red cluster after the next peak. To the topmost right are a number of high reachabilities, which would stand for noise. Note that some of the noise points are classified to be withing a cluster as they are in the density region of some of our clusters.
It is of course easy to make the connections with the color coding and interpreting that without the color codes is much harder but we already see a pattern when compared to the 2D plot: We chose a \(k\) that elicited that we have 5 clusters, that the densest cluster is the green one, followed by black, red and then the two worms. So generally we can say that 1) the deeper a valley the clearer the clustering 2) the more valleys we have the clearer the clustering and 3) the higher the peaks are the clearer the clustering.
Compare this to the following situation where we moved the turqoise cluster closer to the blue one together, shifted the red and black clusters closer to the blue one and increased the variance of the red and back clusters.
x2<-x
x2[which(cl==5),1]<-x2[which(cl==5),1]-0.5
x2[which(cl==5),2]<-x2[which(cl==5),2]-2
x2[which(cl==1),1]<-x2[which(cl==1),1]*1.5
x2[which(cl==1),2]<-(x2[which(cl==1),2]+1.2)*1.5
x2[which(cl==2),1]<-(x2[which(cl==2),1]-1.5)*1.5
x2[which(cl==2),2]<-x2[which(cl==2),2]*1.5
plot(x2, col=cl,pch=19,xlim=c(-3,5),ylim=c(-3,6))
This plot looks less clearly structured (clustered) than the one from before which is easiest to be seen without the coloring
par(mfrow=c(1,2))
plot(x,pch=19, xlim=c(-3,5),ylim=c(-3,6))
plot(x2,pch=19, xlim=c(-3,5),ylim=c(-3,6))
par(mfrow=c(1,1))
Let’s look at what OPTICS tells us about the less clustered situation
We still see that OPTICS does a fairly good job in identifying the number of clusters but notice that the valleys are relatively less deep for black and red and that the peaks are much smaller compared to the first.
Now let’s have turqoise, black and red basically merged.
x3<-x
x3[which(cl==5),1]<-x3[which(cl==5),1]-1
x3[which(cl==5),2]<-x3[which(cl==5),2]-5
x3[which(cl==1),1]<-x3[which(cl==1),1]*2
x3[which(cl==1),2]<-(x3[which(cl==1),2]+1.5)*2
x3[which(cl==2),1]<-(x3[which(cl==2),1]-1.8)*2
x3[which(cl==2),2]<-x3[which(cl==2),2]*2
plot(x3, col=cl,pch=19,xlim=c(-3,5),ylim=c(-3,6))
This plot looks even much less clearly structured (clustered) than the ones from before. This is now just a weird-shaped blob and looks like there’s not really any clustering going on. Again this is easiest to be seen without the coloring.
par(mfrow=c(1,3))
plot(x,pch=19, xlim=c(-3,5),ylim=c(-3,6))
plot(x2,pch=19, xlim=c(-3,5),ylim=c(-3,6))
plot(x3,pch=19, xlim=c(-3,5),ylim=c(-3,6))
par(mfrow=c(1,1))
Let’s look at what OPTICS tells us about the least clustered situation
It can now only identify the dense green cluster and the rest is basically just a blurred mass. Note that the peaks and the valleys are not very different or that the difference of subsequent reachabilities over the ordering is not very large (apart for the green cluster).
Let’s look at the three reachability plots next to each other (note that the highest reachability is largest for the most unclustered situation here, which will be a consideration later when we talk about robustness)
par(mfrow=c(1,3))
plot(ores,col=cl[ores$order],ylim=c(0,2))
plot(ores2,col=cl[ores2$order],ylim=c(0,2))
plot(ores3,col=cl[ores3$order],ylim=c(0,2))
par(mfrow=c(1,1))
So, if we agreee that the first situation is the most clustered, then the second and that one more than the third, we see that in OPTICS this translates to less “up and down” of the reachability plot for the less clustered situations, or less “raggedness” of the reachability plot.
The OPTICS Cordillera
This observation is now at the starting point for the OPTICS Cordillera. Say, we don’t want to look at a reachability plot every time, but want one number that tells us how clustered a plot/data matrix is. We saw that in a more clustered situation, the “up and down” of the reachability plot in the sense of higher peaks and deeper valleys and bigger differences between the two, and the more clusters we have (so more peaks and more valleys), would give us a indication of the clusteredness in the data matrix.
What the OPTICS cordillera now does is precisely that: It measures the raggedness of the reachability plot by taking the norm of the subsequent reachability differences over the OPTICS ordering. The larger norm is, the more peaks, the more valleys we have, the deeper the valleys, the higher the peaks are, in short the more ragged the reachability plot is. And more of that means more clusteredness. This is where the name “cordillera” comes from.
Let’s calculate the OPTICS Cordillera for these situations with the minpts=20, epsilon=10 as before and a maximum winsorization distance of dmax=1.5. This last setting is for robustness so reachabilities larger than 1.5 get set to 1.5 (when looking at the reachability plots above, we see that the highest reachabilities are larger for the situations with less clusteredness because the cluster is much more spread out, see below for more). When comparing the OC over different \(X\), dmax should be fixed to be the same for the OC to be comparable.
cres1<-cordillera(x,minpts=22,epsilon=10,dmax=1.5,scale=FALSE)
summary(cres1)
#>
#> OPTICS Cordillera values with minpts= 22 and epsilon= 10
#>
#> Raw OC: 2.232495
#> Normalization: 83.25
#> Normed OC: 0.24468
cres2<-cordillera(x2,minpts=22,epsilon=10,dmax=1.5,scale=FALSE)
summary(cres2)
#>
#> OPTICS Cordillera values with minpts= 22 and epsilon= 10
#>
#> Raw OC: 1.484021
#> Normalization: 83.25
#> Normed OC: 0.1626477
cres3<-cordillera(x3,minpts=22,epsilon=10,dmax=1.5,scale=FALSE)
summary(cres3)
#>
#> OPTICS Cordillera values with minpts= 22 and epsilon= 10
#>
#> Raw OC: 1.339798
#> Normalization: 83.25
#> Normed OC: 0.146841
We see that both the Raw OC (the length of the “up and down”) and the normalized version normed OC (see below) are highest for the first situation, then for the second, then for the third - exactly as we saw it in the plots above. So the OC is doing what it should.
The OC can also be visualized and is proportional to the fat black enveloping line in the reachability plot (note the effect of dmax which now winsorizes the largest, outlying reachabilities).
We can clearly see that the length of the first black line is the longest (just imagine walking this in the Alps!).
Interpreting and comparing OC values
The raw OC is basically the norm (length) of the up and down between the peaks and the lowest points in each cluster valley. Therefore there is a tendency that the lower minpts is (possibly deeper valleys) and the higher \(n\) is (more possible up and downs) and the higher \(d_{max}\) is (higher possible peaks), the higher the raw OC value becomes. For a single result it is therefore difficult to interpret it in terms of its magnitude due to it capturing a number of aspects such as cluster separation, cluster compactness and number of clusters simultaenously. The normed OC should then be used as it allows for an absolute interpretation as the ratio of how much clusteredness is attained relative to the most clustered arrangement possible given \(k\), \(n\) and \(d_{max}\), or as percent of that if multiplied with 100.
For comparison of different representations of the same data or different data sets with equal and fixed \(k\), \(n\), \(d_{max}\), the raw OC is not dimensionless and can in principle be interpreted in absolute terms between representations of the same data (see above). However, again due to the OC capturing a number of aspects such as cluster separation, cluster compactness and number of clusters simultaenously, it is difficult to say which of these is now the driving factor for the raw OC and we think it is then best to interpret it as a rank (so more raw OC is more clustered, but not exactly how much more). In such a case the normed OC can be interpreted as absolute in terms of its magnitude as it has the same normalization constant for different representations. This allows to add the interpretation of how much of the maximum possible clusteredness is achieved and therefore the values can be interpreted as differences too, so say one representation of a data set with the same \(n\), \(k\), \(d_{max}\) achieves \(5\%\) more clusteredness compared to another representation of the same data set.
For general situations of representation/data sets with different \(n\), \(k\) or \(d_{max}\) for the OC (e.g., different analyses with different parameters, or different data sets) we need to point out that the OC values are typically not directly comparable. In such a situation shoudl the same \(n\), \(k\) or \(d_{max}\) be used for the OC, however, both the raw (again best as a ranking) and the normed value (in the absolute) can be compared. In practice we mostly want to compare different representations of the same data, so often \(n\) and \(k\) are fixed anyway, which usually means we just need to fix \(d_{max}\) for the comparison to be admissable.
The normed OC is further also comparable as a goodness-of-clusteredness statistic between different \(n\), \(k\) and \(d_{max}\) if we want to relatively compare how much more of the most possible clusteredness given the setup is attained for different representations or data sets. This can mean that a situation with a higher raw OC can have less goodness-of-clusteredness if it is further away from the maximum possible clusteredness achievable, which would mean it is relatively less clustered compared to its possible maximum when compared to another.
To sum up we recommend to use the following OC values for different situations
| Single data set with itself (goodness-of-clusteredness) |
normed OC (absolute) |
| Different representations of the same data with the same \(k\), \(d_{max}\) |
raw OC (as ranking), normed OC (absolute) |
| Different data with the same \(n\), \(k\), \(d_{max}\) |
raw OC (as ranking), normed OC (absolute) |
| Different representations of the same data with different \(n\), \(k\), \(d_{max}\) relative to their individual maximum clusteredness |
normed OC (relative) |
| Different data with different \(n\), \(k\), \(d_{max}\) relative to their individual maximum clusteredness |
normed OC (relative) |
| Different representations/data with different \(n\), \(k\), \(d_{max}\) |
not comparable by OC |
Interpretation examples
Let’s look again at the first example. Here we have a single data set that we want to assess for its clusteredness.
We would use the normed OC. For dmax=1.5, we achieve about \(25\%\) of the maximum clusteredness (which would be if we had \(20\) clusters with 22 points each that exactly coincide and each cluster was at a reachability of 1.5 form the neighbouring cluster).
Let’s now compare two representations of the same data once in the original and once where we press the y column together (rescaled to a 10th; not that this would make any sense)
rep2<-x
rep2[,2]<-rep2[,2]/10
par(mfrow=c(1,2))
plot(x,col=cl,asp=1,ylim=c(-3,6),xlim=c(-3,6))
plot(rep2,col=cl,asp=1,ylim=c(-3,6),xlim=c(-3,6))
cres1a<-cordillera(rep2,minpts=22,epsilon=10,dmax=1.5,scale=FALSE)
summary(cres1a)
#>
#> OPTICS Cordillera values with minpts= 22 and epsilon= 10
#>
#> Raw OC: 1.118171
#> Normalization: 83.25
#> Normed OC: 0.1225508
We see that there is lower clusteredness in the second representations (\(12\%\) of achievable).
Customizing the OPTICS Cordillera
The OPTICS Cordillera can be customized with a number of parameters to allow for explore clusteredness in different settings and under different conditions.
To show these effects, we use the following toy data set with a compact cluster of 4 points where two points are denser (nested cluster of 2 points) (Cluster 1 and Cluster 1a nested within), a less compact cluster with three points (Cluster 2), a third cluster with three points of which 2 points are closer (nested cluster) (Cluster 3 and Cluster 3a nested within) and an outlier.
x<-cbind(c(3.7,-4.1,-4.3,-4.6,-4.3,3.2,6.1,4.3,-0.5,-2,-1.7),c(-7,-0.15,-0.1,1.1,-1.2,0.4,1.2,2.5,4.7,6,6.5))
cl<-c("grey",rep("black",4),rep("red",3),rep("green",3))
row.names(x)<-1:dim(x)[1]
plot(x,type="n")
text(x,col=cl,labels=row.names(x),asp=1)
minpts- The most important parameter is probably
minpts and is a parameter for OPTICS (see there for more). It gives the minimum number of points that must make a cluster. It must be at least 2. This will influence how the clusteredness is assessed, for example if we have many 2 point clusters that are well separated and set minpts=2, OC will be high, but if we set it higher, say minpts=3, clusteredness will assessed to be lower because the 2 point clusters are no longer seen as clusters.
- To illustrate for the above data, we have two three point clusters and one 4 point cluster. So if we say
we get a relatively high OC value, as all the clusters are found and the clusters have at least three objects. Compare this to setting minpts=4
which is much lower because the three point clusters no longer count as their own clusters and the reachabilities get extended between Cluster 2 and Cluster 3 so we don’t really have any peaks. The points from 6 to 11 are basically now seen as a single cluster with very low density of at least four points; the only cluster really identified is the 4 point cluster (Cluster 1) with observations 4, 5, 3, 2 that have low density. If we now set minpts to 2, we also allow the OC to take the nested clusters 3a and 1a into account in their own right which prior have been just seen as parts of the larger clusters, which increases the raw OC (but note that the normalization constant is now much higher—in accordance with the explanation in the interpretation section— making the normed OC for \(k=2\) just a little bit higher than in the minpts=3 situation).
epsilon- Is another OPTICS parameter (see there for more) and governs up to which distance we consider points to be neighbours in a cluster. This can be used in two ways: Either to make the OPTICS runtime quicker, which can be beneficial with a lot of points and, more importantly, to designate points as noise if the density in their vicinity is not large enough (in the OC, a noise point gets assigned \(d_{max}\) and they are ordered next to each other so their reachbility difference is zero). In our situation we have a point that is dedicated as noise (point 1) and in the setup before (at least three points per cluster) it would therefore be included as a possible point in cluster 2 where it is in the neighbourhood (density-reachable) of the observations 6, 7, 8. This would however make that cluster very non-dense, giving high reachability from 7 (the succesor in the ordering) to that point and overall decreasing the OC.
We might not want that because we know it’s an outlier or that there is noise in our data, so we reduce epsilon to 7 and the point no longer becomes “density-reachable” from the other points. Now the point is no longer in any cluster and this makes the individual clusters comparatively more dense, meaning the OC increases
Typically, though, this is a fickle parameter and it might be best to just leave it large enough to possibly include all points (that was OPTICS improvement over DBSCAN, to abstract from having to set \(\varepsilon\). We also caution against setting epsilon too low, which can assign too many points as noise (optics reachability of \(\infty\), or OC reachability of \(d_{max}\)), e.g.
dmax and rang- These are parameters that allow for robustness in the OC. In the Cordillera, all reachabilities that are larger than
dmax get winsorized to dmax (including points with a reachability of inf in OPTICS). For the normalization, dmax also plays a role as it is the maximum reference distance that needs to exist between neighbouring clusters in a perfect clusteredness situation. dmax should be larger than most reachabilities but smaller than any outlying reachability to not inflate the OC. For example in the situation with low clusteredness above (x3)
We have a relatively large reachability (around 3) for the last point in the ordering. But the distribution of reachabilities is highly right-skewed and \(90\%\) of the reachabilities are smaller than 0.54, when we look at the quantiles of the reachability distribution.
: The OC is by itself not robust to such outliers, so with the winsorization it allows to exclude the effect of these outliers. The default dmax is third quantile + 1.5 times the interquartile range, so on in our case here
this would be 0.64. Let’s compare the OC in these two cases
: where we see that the raw OC is much larger in the first case due to the outlier. This is a bit mitigated in the normed OC, but we still think that the robust version gives a more useful normalized version too. How to set dmax is of course up to the user, but it shouldn’t be too low either (e.g. set at the median because) this can result in too low a value (or even 0). The rang is the interval from any lower bound to any upper bound of the reachability and allows to fine tune the relative contribution for the OC if we want to measure the ups and downs not against zero, but against any other lower/upper bound for the reachabilities (e.g., the minimum and maximum reachability). We don’t expect users to use this parameter, though, and it’s default is (0,dmax). If it is given however and the maximum doesn’t agree with dmax, dmax is overridden, so use this with extreme caution.
q- Is a parameter relating to the OC and gives the \(L_p\)-space for the norm that is used for the OC. The most important values are
q=1 where the sum of absolute reachability differences are taken over the ordering and q=2 where the Euclidean norm is used, so the square root of the sum of squared reachability differences over the ordering. For consistency, we believe q=1 is a good idea if the reachabilities are calculated from a matrix of Manhattan distances q=2 (the default) when calculated from a Euclidean distance matrix.
distmeth- Allows to calculate a distance matrix for the \(X\) argument if it is not already a distance matrix. It supports all distances of
stats::dist. Default is “euclidean”.
scale- this allows to scale the \(X\) argument before calculating the cordillera.
scale=FALSE or scale=0 takes \(X\) as it is. If scale=TRUE or 1, standardisation is to mean=0 and sd=1 for each columns. If scale=2, no centering is applied and scaling of each column is done with the root mean square of each column. If scale=3, no centering is applied and scaling of all columns is done as X/max(standard deviation(allcolumns)). If scale=4, no centering is applied and scaling of all columns is done as X/max(rmsq(allcolumns)). Default is scale=FALSE. The scaling options are typically best used if we want to compare representations of the same data that can have diffferent scales and where we are only interested in relative distances between points, as in most dimensions reduction methods.