# Inferring hidden categorical variables from multimodal response data

Lets say I have a response variable that is bimodally distributed. Below is a simple example with 1 predictor variable:

This can be dealt with by adding a categorical variable that splits the data into two groups. For example, in the toy data above you could imagine y=child’s height (full grown) and x=parent’s height. The data above the best fit line would be males and the data below the best fit line would be females. Let’s say, for arguments sake, that we didn’t originally record the gender of each child.

First (small) question: I suppose one could use a clustering algorithm to recover the two groups in the data, and then create a linear regression model with the clustering output as an additional categorical predictor variable. Is this ever done? Are there problems with this logic?

Second (more important) question: Let’s say we want to predict the height of a new child, but we don’t know his/her gender. Rather than giving the expected value (red line), I want to return two possible values (i.e. height if male, or height if female). What is the best way to do this?

Preliminary answer: Use the clustering approach I mentioned above, and fit an expanded linear regression model with an additional categorical variable which nominally refers to gender. Then, given a new observation x, predict y in the expanded model for both possible cases.

Difficulties: I am interested in how to solve this problem for much higher dimensional data sets, and specifically how to automate it. In general, it will not be known how many clusters there are in the data, or which variables to cluster the data along. The data I am interested in might also be highly non-linear. Consider the example below:

I think this is an interesting question for dynamical systems that switch between multiple attractor states. For example a neuronal network with a noisy (and uncontrolled) input that can “bump” or switch the activity pattern between multiple (unknown) states. This might lead to different clusters in the response variable as in the above example. See: http://www.scholarpedia.org/article/Attractor_network

Matlab code for generating figure 2:

N = 2e2;
x = rand(N,1);
y1 = nan(N,1); y2 = nan(N,1);
y1(x<0.4) = 1.2*x(x<0.4)+0.2+0.03*randn(sum(x<0.4),1);
y2(x<0.4) = x(x<0.4)+0.6+0.03*randn(sum(x<0.4),1);
y1(x>=0.4 & x<0.6) = 0.5*x(x>=0.4 & x<0.6)+0.01*randn(sum(x>=0.4 & x<0.6),1);
y2(x>=0.4 & x<0.6) = 0.6*x(x>=0.4 & x<0.6)+0.01*randn(sum(x>=0.4 & x<0.6),1);
y1(x>=0.6) = 0.5+0.02*randn(sum(x>=0.6),1);
y2(x>=0.6) = 0.5+0.02*randn(sum(x>=0.6),1);
y = [y1; y2];
x = repmat(x,2,1)+0.02*randn(2*N,1);
plot(x,y,'ob')
xlabel('x','FontSize',20)
ylabel('y','FontSize',20)

• You are looking for model-based clustering techniques, such as finite mixture models. Dec 15 '13 at 19:37
• There is nice DP-means clustering code for R that can be found here: johnmyleswhite.com/notebook/2012/06/26/… Dec 15 '13 at 20:04
• @AlexWilliams if you share the data used for those plots I can try to answer your questions directly using the method in my answer below. Dec 16 '13 at 17:10
• @Flask -- Thanks for the response. I will post the code soon. After you generate the clusters - would you then fit a separate statistical model to each one? Dec 16 '13 at 19:25

I made a 3d version of John Myles White's nice code he shared here: http://www.johnmyleswhite.com/notebook/2012/06/26/bayesian-nonparametrics-in-r/

The only difference is the extra dimension and I added a normalization step for each dimension. This algorithm is nice because rather than setting the number of clusters beforehand you choose a parameter "lambda" that determines the minimum distance between a point and the center of a cluster. Any points farther than this distance will start new clusters. I see no reason you couldn't simply add additional dimensions as needed although it is somewhat tedious with this code.

Once you have assigned a cluster to each data point you can fit lines/whatever to each separately.

Modified data generating function:

generate.data <- function()
{
n <- 100
data <- data.frame(x = rep(NA, n), y= rep(NA, n), z= rep(NA, n))
mu.x <- c(0, 5, 10, 15)
mu.y <- c(0, 5, 0, -5)
mu.z <- c(0, 5, 0, -5)

for (i in 1:n)
{
assignment <- sample(1:4, 1)
data[i, 'x'] <- rnorm(1, mu.x[assignment], 1)
data[i, 'y'] <- rnorm(1, mu.y[assignment], 1)
data[i, 'z'] <- rnorm(1, mu.z[assignment], 1)
}

return(data)
}


Modified clustering algorithm:

dp.means.3d <- function(data, lambda = 1, max.iterations = 100, tolerance = 10e-3){

data$x<-(data$x-mean(data$x))/(sd(data$x))
data$y<-(data$y-mean(data$y))/(sd(data$y))
data$z<-(data$z-mean(data$z))/(sd(data$z))

n <- nrow(data)
k <- 1
assignments <- rep(1, n)
mu.x <- mean(data$x) mu.y <- mean(data$y)
mu.z <- mean(data$z) converged <- FALSE iteration <- 0 ss.old <- Inf ss.new <- Inf while (!converged && iteration < max.iterations) { iteration <- iteration + 1 for (i in 1:n) { distances <- rep(NA, k) for (j in 1:k) { distances[j] <- (data[i, 'x'] - mu.x[j])^2 + (data[i, 'y'] - mu.y[j])^2 + (data[i, 'z'] - mu.z[j])^2 } if (min(distances) > lambda) { k <- k + 1 assignments[i] <- k mu.x[k] <- data[i, 'x'] mu.y[k] <- data[i, 'y'] mu.z[k] <- data[i, 'z'] } else { assignments[i] <- which(distances == min(distances)) } } for (j in 1:k) { if (length(assignments == j) > 0) { mu.x[j] <- mean(data[assignments == j, 'x']) mu.y[j] <- mean(data[assignments == j, 'y']) mu.z[j] <- mean(data[assignments == j, 'z']) } } ss.new <- 0 for (i in 1:n) { ss.new <- ss.new + (data[i, 'x'] - mu.x[assignments[i]])^2 + (data[i, 'y'] - mu.y[assignments[i]])^2 + (data[i, 'z'] - mu.z[assignments[i]])^2 } ss.change <- ss.old - ss.new ss.old <- ss.new if (!is.nan(ss.change) && ss.change < tolerance) { converged <- TRUE } } centers <- data.frame(x = mu.x, y = mu.y, z = mu.z) return(list("centers" = centers, "assignments" = factor(assignments), "k" = k, "iterations" = iteration)) }  Generate Data and find the clusters: #install.packages("rgl") require(rgl) dat<-generate.data() clusters<-dp.means.3d(dat) plot3d(dat, col=clusters$assignments, size=10)


Edit:

The Results:

So, as you can see it doesn't like that the results are not multivariate normal for larger lambda.

Here is the part of the code that determines when to make a new cluster:

  for (j in 1:k)
{
distances[j] <- (data[i, 'x'] - mu.x[j])^2 + (data[i, 'y'] - mu.y[j])^2
}


There may be a better way if we want to make no assumptions about the data at all. Perhaps it should also include a step that attempts to minimize the maximum distance to other points part of the same cluster.

> fits
Lambda Cluster   n  Intercept        Slope SumSquared       RMSE
[1,]      1       1  98 0.47108474  0.037932047 0.03386304 0.01858874
[2,]      1       2  32 0.28115152  0.803120715 0.02721487 0.02916273
[3,]      1       3  66 0.01454559  0.523555233 0.06129777 0.03047547
[4,]      1       4  76 0.50312253 -0.004766139 0.02985764 0.01982079
[5,]      1       5  29 0.23240685  0.750050744 0.02459047 0.02911953
[6,]      1       6  19 0.30395055  0.944757052 0.01314973 0.02630762
[7,]      1       7  32 0.71040200  0.446582407 0.04330122 0.03678537
[8,]      1       8  27 0.61745795  0.637197982 0.02369394 0.02962353
[9,]      1       9  21 0.77826596  0.493610959 0.01757857 0.02893224
[10,]      2       1  37 0.24003453  1.101926169 0.03441819 0.03049953
[11,]      2       2  66 0.01454559  0.523555233 0.06129777 0.03047547
[12,]      2       3 174 0.49416121  0.005600235 0.06415743 0.01920211
[13,]      2       4  43 0.22599167  0.919170075 0.04074510 0.03078247
[14,]      2       5  41 0.68403859  0.734703194 0.05907046 0.03795712
[15,]      2       6  39 0.61635785  0.736739137 0.04078044 0.03233655
[16,]      3       1 174 0.49416121  0.005600235 0.06415743 0.01920211
[17,]      3       2 127 0.36147381 -0.130376879 0.94977762 0.08647869
[18,]      3       3  99 0.65870905  0.449620540 1.29255466 0.11426333
[19,]      4       1 174 0.49416121  0.005600235 0.06415743 0.01920211
[20,]      4       2 127 0.36147381 -0.130376879 0.94977762 0.08647869
[21,]      4       3  99 0.65870905  0.449620540 1.29255466 0.11426333
[22,]      5       1 174 0.49416121  0.005600235 0.06415743 0.01920211
[23,]      5       2 127 0.36147381 -0.130376879 0.94977762 0.08647869
[24,]      5       3  99 0.65870905  0.449620540 1.29255466 0.11426333
[25,]      6       1 300 0.27214072  0.247572658 1.96115772 0.08085291
[26,]      6       2 100 0.65807225  0.441094682 1.35260593 0.11630159


DP-means 2d algo:

dp.means <- function(data, lambda = 1, max.iterations = 100, tolerance = 10e-3)
{

data$x<-(data$x-mean(data$x))/sd(data$x)
data$y<-(data$y-mean(data$y))/sd(data$y)

n <- nrow(data)
k <- 1
assignments <- rep(1, n)
mu.x <- mean(data$x) mu.y <- mean(data$y)

converged <- FALSE
iteration <- 0

ss.old <- Inf
ss.new <- Inf

while (!converged && iteration < max.iterations)
{
iteration <- iteration + 1

for (i in 1:n)
{
distances <- rep(NA, k)

for (j in 1:k)
{
distances[j] <- (data[i, 'x'] - mu.x[j])^2 + (data[i, 'y'] - mu.y[j])^2
}

if (min(distances) > lambda)
{
k <- k + 1
assignments[i] <- k
mu.x[k] <- data[i, 'x']
mu.y[k] <- data[i, 'y']
} else
{
assignments[i] <- which(distances == min(distances))
}
}

for (j in 1:k)
{
if (length(assignments == j) > 0)
{
mu.x[j] <- mean(data[assignments == j, 'x'])
mu.y[j] <- mean(data[assignments == j, 'y'])
}
}

ss.new <- 0

for (i in 1:n)
{
ss.new <- ss.new + (data[i, 'x'] - mu.x[assignments[i]])^2 + (data[i, 'y'] - mu.y[assignments[i]])^2
}

ss.change <- ss.old - ss.new
ss.old <- ss.new

if (!is.nan(ss.change) && ss.change < tolerance)
{
converged <- TRUE
}
}

centers <- data.frame(x = mu.x, y = mu.y)

return(list("centers" = centers, "assignments" = factor(assignments), "k" = k, "iterations" = iteration))
}


Generate Data:

N = 2e2
x = runif(N,0,1)
y1 = matrix(nrow=N,ncol=1)
y2 = matrix(nrow=N,ncol=1)
n1<-length(which(x<0.4))
y1[which(x<0.4)] = 1.2*x[which(x<0.4)]+0.2+0.03*rnorm(n1,0,1)
y2[which(x<0.4)] = x[which(x<0.4)]+0.6+0.03*rnorm(n1,0,1)

n2<-length(which(x>=0.4 & x<0.6))
y1[which(x>=0.4 & x<0.6)] = 0.5*x[which(x>=0.4 & x<0.6)]+0.01*rnorm(n2,0,1)
y2[which(x>=0.4 & x<0.6)] = 0.6*x[which(x>=0.4 & x<0.6)]+0.01*rnorm(n2,0,1)

n3<-length(which(x>=0.6))
y1[which(x>=0.6)] = 0.5+0.02*rnorm(n3,0,1)
y2[which(x>=0.6)] = 0.5+0.02*rnorm(n3,0,1)
y = c(y1, y2)
x = rep(x,2)+0.02*rnorm(2*N,0,1)

dat<-data.frame(cbind(x,y))


Find Clusters and fit lines:

fits=NULL
par(mfrow=c(3,2))
for(lambda in 1:6){
clusters<-dp.means(dat, lambda=lambda, max.iterations = 10000, tolerance = 10e-3)

plot(x,y, xlab="x", ylab="y", pch=16,
col=rainbow(clusters$k)[clusters$assignment],
main=paste("Lambda=",lambda))
legend("topright", legend=paste("Cluster",1:clusters$k), col=rainbow(clusters$k), pch=16, pt.cex=2, ncol=2)

for(i in 1:clusters$k){ sub<-dat[which(clusters$assignments==i),]
fit<-lm(sub$y~sub$x)
n<-length(sub$x) SS<-sum((fit$residuals)^2)
rmse<-sqrt(SS/n)
lines(sub$x,fit$fitted.values, lwd=2)
fits<-rbind(fits,cbind(lambda,i,n, rbind(fit\$coefficients), SS, rmse))
}

}
colnames(fits)<-c("Lambda","Cluster","n","Intercept","Slope","SumSquared","RMSE")

}