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

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?

• I'd say they form much stronger vertical than horizontal lines. SCNR. And can't you just use the number of values you have for each x value as the number of clusters? – Has QUIT--Anony-Mousse Jul 6 '12 at 16:10
• The problem is not sufficiently well defined yet. Additional information that would be useful can include (1) what might be known about how and why the coordinates vary; (2) whether "horizontal" can include slopes that deviate slightly from horizontal; (3) the purpose of this analysis; (4) typical sizes of datasets that will be encountered; and (5) any special restrictions, such as whether there will be the same number of points associated with each $x$ value. – whuber Jul 6 '12 at 16:14
• @whuber: 1. the variation in y is due to random noise. 2. no. 4. very large, say 200 lines with 1000 points each. 5. usually yes, but some x values may have a few (say 1->10) y values missing. – Andy Jul 6 '12 at 16:35
• In that case it sounds like you could select the minimum value for each x, remove any high outliers relative to those data, and call what you have a horizontal line. Repeat until done. – whuber Jul 6 '12 at 16:49
• re missing values: The missing values cause high outliers to appear previously, which you exclude. re geophysics: sedimentary strata are almost never perfectly horizontal (or uniformly thick)! You may be making too strong an assumption. – whuber Jul 6 '12 at 17:00

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.

• Thanks Anony-Mousse! Starting from the column with the lowest number of peaks sounds like a good idea since this will increase the chanches of the neighbour columns having a y value on this line. But still this may happen in which case this column must be skipped. Minimizing the cost function sounds like a good idea, but when growing from one column there may be a connection between two other columns that have a lower cost. Would it not be better (if it is possible) at each step to join the two neighbouring columns with the lowest cost? – Andy Jul 6 '12 at 17:36
• In which case the algorithm starts to sound a lot like hierarchical agglomerate clustering. – Andy Jul 6 '12 at 17:37
• BTW. The cost between connected points may also better than the cost between each point and the mean (k-means) since this would relax the requirement that the lines are horizontal (?). Only thing is that I want minimize this globally: summed over all line, summed over all connections. Therefore I think that a global greedy algorithm may give better results than a local one? – Andy Jul 6 '12 at 17:48
• Well, but doesn't the x axis indicate a spatial relationship? That is why I would progress along x. The final optimization passes can still correct errors that you made in the beginning, too! – Has QUIT--Anony-Mousse Jul 6 '12 at 18:26
• Global may yield a numerically lower score, but it makes the result much less meaningful. The key to a good result is to put the knowledge that you have into model constraints! – Has QUIT--Anony-Mousse Jul 6 '12 at 18:28

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.

• Thanks whuber, the results looks very impressive. I am not familiar with Lowess smoothing. From what you mention it sounds like it helps both with finding a (segmentwise?) straight line and with identifying outliers? If so over how large a window do you check for outliers? Because it is some form of local regression? – Andy Jul 7 '12 at 17:13
• Andy: Yes, it is locally weighted regression. I found that for these simulated data a window of about 20% of all the x-values gave good results. But just about any robust local smoother ought to work, such as Tukey's 3RSSH. – whuber Jul 9 '12 at 12:18