How do I cluster data where some observation can not be in the same cluster?

Question

Say I have 12 (x,y) positions:

data =
    1.0000    0.8500
    2.0000    1.1000
    3.0000    1.2000
    4.0000    0.9000
    1.0000    2.0000
    2.0000    2.1000
    3.0000    1.9000
    4.0000    2.0000
    1.0000    3.1000
    2.0000    3.1000
    3.0000    2.9000
    4.0000    2.8000

I want to classify these points into horizontal lines. From a scatterplot it is easy for a human to see that the points form 3 horizontal lines; one around y=1, another around y=2 and the last around y=3:

Which classification algorithm should I use to classify the points?

The easiest algorithm seem to be k-means. One problem with this is that I need to specify the number of clusters. How do one typically get around this problem? I was thinking that I may run it with an increasing number of clusters, n, and calculate the total square deviation from cluster mean and once the decrease in total square deviation relatively to the last number of clusters, n-1, is less than a certain threshold (say 1%) then n is a sufficient number of clusters.

Another algorithm I was considering using is hierarchical agglomerate clustering. The solution here would form a dendrogram, and like for the k-mean I may choose a level with n clusters for which the total deviation is below some threshold criteria.

One important piece of information that I do not understand how to use though in either of these algorithms is this; no cluster may have the same x position more than once! So far I have just ignored the x coordinate and clustered based on y values.

How can I use the criteria that no cluster may have the same x position more than once with a cluster algorithm?

Thanks in advance for any answers!

Answer 1

I assume the x values are somewhat controlled and not really measurements. You shouldn't assume these to be "equivalent" to the y values, so stay away from k-means, which assumes that Euclidean distance is a sensible distance measurement.

Instead, treat the x axis as a source of given similarity and the y axis as a kind of measured similarity. Try to see this as time series.

The most simple method would be just to always align one "column" with the next. For example by sorting them by y value, and mapping them in an ordered way. Lowest to lowest, second lowest to second lowest etc. Maybe introduce some threshold.

Do not just throw the data into some random clustering algorithm, this will not get you anywhere. Try to find a good model for your data, and try to optimize this.

A more advanced approach would be to:

• Find the column with the lowest number of peaks.
• Grow both to the left and right, column by column
• For matching the next column, start with matching the median to the median, then align them one by one. Define the cost function $$C=\sum_i |y_{old}(i) - y_{new}(i)|$$ Try to optimize this cost function by moving the alignments one step at a time, but never criss-crossing lines and maybe also never losing lines!
• Proceed to the next columns outward.

At the very end, maybe try another pass at optimizing the cost function outward-in. And add partial lines that disappeared somewhere.

Note that on your example data set, the result will be trivial: it will map the lowest value for each x to the lowest of each other x, the highest to the highest each, and so on. It's a way too simple example.

Answer 2

One approach is to identify strata from the bottom up (or top down), "peeling away" the data for each stratum as one proceeds.

Given the nature of geological data--strata are almost never perfectly horizontal, they vary in thickness, and so on--it seems best to use a robust non-parametric method to fit the elevations locally. Lowess will do this pretty well.

As an example, I generated 20 layers of 100 points at uniformly random distances among the layers and a uniform underlying slope of -1:100, applying some normally distributed variation to their elevations. I then eliminated 60 of those points at random (an average of three per layer).

(The close spacings in the 5-6 range of elevations will confound any standard clustering procedure I know of: this is intended to be a challenging dataset.)

After using Lowess to smooth each layer, any high outliers were removed and put back into the remaining dataset for future processing. This is not ideal: instead, this step should be provisional. After processing the next layer, such outliers ought to be re-examined to assess which of the two layers it more properly belongs to: the current one, the next one, or (possibly) an even higher one.

(Outliers were identified using a simple robust technique that sets an upper "fence" based on the median and an upper percentile of the Lowess residuals, akin to how an upper fence is established for a boxplot.)

Colors reflect the classifications of the data. Lines show the Lowess smooth of each set of classified points. The number of groups ("clusters") was automatically determined by the algorithm, which stops when it runs out of points.

A close look shows a few misclassifications due to this non-ideal approach. Even so, the contacts themselves (shown by the Lowess lines) are robust to these mis-classifications; one could post-process the results to correct almost all the errors. The algorithm appears to have done a remarkably good job with closely adjacent contacts (look around heights 5-6, for example). In particular, the smooths never cross (even though there is nothing to prevent this), which is essential for geological meaningfulness.

Because Lowess works in multiple dimensions, the same algorithm would work for spatially extensive data (instead of data along a cross-section as shown here).

The use of Lowess is not essential: any strong, robust, local smoother should work well. Indeed, I would have preferred iteratively reweighted least squares, but a good implementation was not ready at hand for making this illustration.