3
$\begingroup$

Today I faced what I think is a very simple problem, but could't solve. I have this plot (data is below)

with(mydata, plot(x, y))

enter image description here

It's clear that there are two groups here, right? I would like to find a line that discriminates these two groups, like this one (found manually).

plot(mydata$x, mydata$y); abline(9, .36)

enter image description here

Questions

What technique should I use? Could anyone provide R code to solve it, please?

mydata <- structure(list(y = c(24.91543, 20.93154, 33.38601, 32.6498, 38.96774, 20.51704, 10.64043, 35.69466, 28.72396, 31.80371, 34.11269, 31.64434, 25.5895, 33.58181, 34.34093, 25.92858, 32.21135, 30.19319, 30.27322, 24.43336, 26.00838, 35.42962, 26.46969, 32.88677, 21.21332, 25.6222, 36.60242, 27.00271, 40.59243, 29.02613, 30.10023, 15.59976, 41.07338, 31.13055, 34.00474, 37.13078, 30.70407, 30.80497, 32.9493, 42.1743, 38.7733, 26.01051, 32.47489, 21.95968, 25.1595, 38.58872, 27.98961, 42.83003, 31.08639, 34.62195, 38.96774, 32.56765, 16.40093, 35.18166, 39.06299, 34.79893, 22.88687, 35.44114, 37.627, 34.91972, 25.11421, 26.3312, 42.68491, 34.94547, 16.22444, 29.20567, 21.39507, 27.31063, 38.83989, 18.57524, 28.76647, 39.07916, 36.94783, 25.35636, 38.4268, 42.10769, 44.762, 21.71412, 36.63178, 25.50169, 37.27536, 28.62975, 25.89833, 26.34055, 20.88283, 29.91666, 40.16116, 38.72035, 30.23885, 29.58979, 19.41681, 23.4266, 31.84684, 32.46026, 34.90148, 23.12917, 30.60789, 15.79587, 33.64789, 30.11133), x = c(37.139999, 50, 60, 90, 96.669998, 46.669998, 20, 62.860001, 37.139999, 54.290001, 93.330002, 48.57, 70, 90, 60, 73.330002, 54.290001, 45.709999, 40, 28.57, 70, 60, 34.290001, 51.43, 50, 37.139999, 68.57, 37.139999, 80, 34.290001, 37.139999, 30, 80, 83.330002, 54.290001, 65.709999, 42.860001, 83.330002, 48.57, 85.709999, 74.290001, 70, 54.290001, 22.860001, 31.43, 68.57, 76.669998, 82.860001, 42.860001, 57.139999, 96.669998, 51.43, 26.67, 60, 96.669998, 60, 28.57, 62.860001, 65.709999, 62.860001, 28.57, 70, 85.709999, 62.860001, 30, 80, 28.57, 73.330002, 71.43, 50, 48.57, 77.139999, 93.330002, 70, 68.57, 85.709999, 100, 22.860001, 93.330002, 70, 68.57, 80, 66.669998, 70, 14.29, 42.860001, 96.669998, 96.669998, 83.330002, 51.43, 14.29, 60, 45.709999, 51.43, 57.139999, 28.57, 45.709999, 23.33, 54.290001, 83.330002)), .Names = c("y", "x"), row.names = c(NA, -100L), class = "data.frame")

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ What do you mean "find the line?" Do you have the intercept and slope? What you are probably looking for is discriminant analysis. statmethods.net/advstats/discriminant.html $\endgroup$
    – le_andrew
    Jun 19, 2015 at 16:03
  • $\begingroup$ I would like to find the intercept and the slope. I can't use discriminant analysis because I don't have the groups yet, am I right? $\endgroup$
    – Julio Trecenti
    Jun 19, 2015 at 16:05
  • $\begingroup$ Maybe what I want is an unsupervised algorithm that separates space in two groups, like kmeans. However, kmeans does not provide a function, just centroids. I've already ran it but got no luck (it separated my plot in left and right)... $\endgroup$
    – Julio Trecenti
    Jun 19, 2015 at 16:12

3 Answers 3

Reset to default
2
$\begingroup$

Not sure if there is a "better" solution, but you'll get what you want from a mixture of clustering and ANCOVA.

First, create a 2 group cluster (using, for instance, Kmeans). Then use the cluster assignments as a categorical predictor in an ANCOVA. The intercept and the slope for your continuous variable are the values you are looking for (assuming you use sum-to-zero codes on your categorical predictor).

set.seed(1234)

x1<-rnorm(50,0,.1)
shift<-c(rep(0,25),rep(3,25))
y<-2+3*x1+shift+rnorm(50,0,.3)

fit<-kmeans(cbind(scale(x1),scale(y)),2)

x2<-factor(fit$cluster)

contrasts(x2)<-contr.helmert(2)

plot(x1,y,col=fit$cluster)

sum.1<-summary(lm(y~x1+x2))

abline(sum.1$coefficients[1],sum.1$coefficients[2])
sum.1

Admittedly, kmeans clustering may not always give you exactly the results you are looking for. A nearest neighbor hierarchical clustering may be better.

Improve this answer
$\endgroup$
5
  • $\begingroup$ This assumes the slopes of the two lines are the same, possibly an interaction between x1 and x2 would be better? $\endgroup$
    – Matthew Drury
    Jun 19, 2015 at 16:32
  • $\begingroup$ Very nice. I think that the "same slope" assumption works great for my data. It is clear now that I need some usupervised + some regression to solve this kind of problem. Thanks! $\endgroup$
    – Julio Trecenti
    Jun 19, 2015 at 16:37
  • 1
    $\begingroup$ This is where the question starts to get too vague, Matthew. Technically there are infinite lines that separate them, assuming the sets are convex. I think the general case can be solved using separating hyperplanes and support vectors, but its been WAY too long since I've done that kind of math. Even with different slopes, though, the ANCOVA results should find A line that separates them, though maybe its not the "best" however you define that. $\endgroup$
    – le_andrew
    Jun 19, 2015 at 16:52
  • $\begingroup$ E.g. if you do the same procedure with the following y values (with different slopes for the different groups) you'll still get a separating line, even though it looks odd. y2<-2+3*x1*shift+shift+rnorm(50,0,.3) $\endgroup$
    – le_andrew
    Jun 19, 2015 at 16:54
  • $\begingroup$ As Vincent demonstrates in his answer below (and you mention in your answer), simple kmeans would perform quite poorly here (due to isotropic covariance). Mixture of Gaussians should work well. I don't think hierarchical clustering, just based on Euclidian distances, would work either. $\endgroup$
    – Zhubarb
    Jun 22, 2015 at 12:21
8
$\begingroup$

In order to find this line, you have first to find the two groups. There is a super cool R library called FlexMix that allows you to find the two clusters. The method used here aims at finding 2 (k=2) different linear models (y~x) in the data. The methods is based on an iterative process that can yield slightly different results between 2 different runs, that's why I used set.seed(123).

set.seed(123)
require(flexmix)
res.flex <- flexmix(y~x, data=mydata, k=2)  
plot(y~x, data=mydata, col=res.flex@cluster)

Flexmix Result

Once you have the groups, you can apply your favorite classification algorithm, like discriminant analysis:

require(MASS)
Z <- scale(mydata)
g <- factor(res.flex@cluster)
res.lda <- lda(g~., data=data.frame(g=g, Z))

gmean <- res.lda$prior%*%res.lda$means
const <- drop(gmean%*%res.lda$scaling)

slope <- -res.lda$scaling[1]/res.lda$scaling[2]
intercept <- const/res.lda$scaling[2]

The resulting linear separator is in green in the following plot.

plot(y~x, data=Z, col=g)
abline(intercept, slope, col=3)

Linear Separator

Just for fun, I tried to find the two groups with K-Means, and the resulting groups was not great at all :

Kmeans

EDIT 1: Following @le_andrew's comment below, here are the resulting clusters after scaling:

Z <- scale(mydata)
plot(y~x, data=Z, col=kmeans(Z,2)$cluster)

After Scaling

EDIT 2: I found in these examples how to compute the slope and intercept of the linear separator, I directly edited my previous code.

EDIT 3: Since flexmix uses posterior probabilities to identify the clusters, maybe we could use its output to directly deduce the linear separator. I don't know yet how to do that, though.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Since the scales of your data are different, you should be standardizing before doing kmeans. $\endgroup$
    – le_andrew
    Jun 19, 2015 at 17:08
  • 2
    $\begingroup$ Without more guidance to the algorithm, who's to say that the FlexMix grouping is "better" than the K-Means (with or without scaling)? Different, yes, but better? Better because that is the grouping the OP wants? If the OP already knows what he wants, he could just draw the line himeself, as he did. $\endgroup$
    – Mark L. Stone
    Jun 22, 2015 at 13:15
  • 1
    $\begingroup$ This is a very good remark. FlexMix was developed to identify a mixture of models. I used it here because there seems to be a mixture of two linear models in the scatterplot and because none of the methods I tried before using it could recover the two groups the OP wanted to automatically get. If you look closely at the OP's Figure and the result of FlexMix, you will see that the two clusters are slightly different. $\endgroup$
    – Vincent Guillemot
    Jun 22, 2015 at 13:29
0
$\begingroup$

You could call the linear model and get the intercept and slope from the summary, as the coefficients give you the slope. 

lm(y ~ x)

This should give you the coefficients (i.e. the slope).

Improve this answer
$\endgroup$
1
  • $\begingroup$ It's just a regression line, it does not separate my plot in the "right" two groups. $\endgroup$
    – Julio Trecenti
    Jun 19, 2015 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.