Very steep decrease in information criteria for mixture models with more components

Question

I am analyzing data using mixture modeling. When I plot the information criteria (the BIC) for a series of models (with different model specifications and different number of mixture components), I notice a very steep decrease in the BIC. E.g., in the photo below, the fitted model represented by the grey lines (and points) with five profiles (or components) has a much lower BIC than that with four components (the BIC decreases by around 50%). Similar for the model represented by the red line between five and six mixture components.

Is this a cause for concern? Or is it likely that this is due to real improvements in the model fit when an additional component is added? It seems like a change of this magnitude, all other things being equal, is cause for concern. What would you to do to follow-up on this analysis to determine whether the change in the log-likelihood is spurious and problematic or indicative of (practically) meaningfully better fitted models?

Answer 1

Here is how to get fractional membership.

library(MASS)

set.seed(1)

mu1 <- c(0,0)
sig1 <- matrix(data=c(1,0,0,1), nrow=2,ncol=2)

mu2 <- c(3,4)
sig2 <- matrix(data=c(1,0,0,1), nrow=2,ncol=2)

mu3 <- c(9,4)
sig3 <- matrix(data=c(1,0.1,0,1), nrow=2,ncol=2)

mu4 <- c(12,0)
sig4 <- matrix(data=c(1,0,0,1), nrow=2,ncol=2)

mu5 <- c(6,-3)
sig5 <- matrix(data=c(1,0.1,0.1,1), nrow=2,ncol=2)

x1 <- mvrnorm(n = 1000, mu = mu1, Sigma = sig1)
x2 <- mvrnorm(n = 1000, mu = mu2, Sigma = sig2)
x3 <- mvrnorm(n = 1000, mu = mu3, Sigma = sig3)
x4 <- mvrnorm(n = 1000, mu = mu4, Sigma = sig4)
x5 <- mvrnorm(n = 1000, mu = mu5, Sigma = sig5)


x <- rbind(x1, x2, x3, x4, x5)
y <- rbind(rep.int(1,times = 1000),
           rep.int(2,times = 1000),
           rep.int(3,times = 1000),
           rep.int(4,times = 1000),
           rep.int(5,times = 1000))

c2 <- kmeans(x=x,centers = 2)

plot(x)
grid()
idx <- which(c2$cluster==1)
points(x=x[idx,1],y=x[idx,2], pch=19, col="Green")
points(x=x[-idx,1],y=x[-idx,2], pch=19, col="Red")

Yielding:

So here is code to estimate AIC (a cousin of BIC):

kmeansAIC = function(fit){

     m = ncol(fit$centers)
     n = length(fit$cluster)
     k = nrow(fit$centers)
     D = fit$tot.withinss
     return(D + 2*m*k)
}


my_bic <- numeric(length = 10)
for(i in 1:10){
     c2 <- kmeans(x=x,centers = i)
     my_bic[i] <- kmeansAIC(c2)

}

plot(log10(my_bic))
points(5,log10(my_bic[5]), pch=19, col="Red")

and here is the result:

This is why I like the "scree plot" approach. I know it is 5 clusters. You know it is 5 clusters. Looking at the differences, the 5-component model is indicated.