Mixture of beta distributions: full example

Question

I have a vector of measured proportions and want to

1. test if these proportions follow a bimodal distribution
2. characterize the two underlying distributions
3. determine for each data point which distribution it belongs to
4. overaly the density distributions on the histogram

I've been using "betamix" from the betareg package but haven't figured out steps 2 & 3. I've searched stackoverflow for solutions but haven't found a clear answer.

Here's my current code:

# parameters of distribution #1 
alpha1 <- 10
beta1 <- 30

# parameters of distribution #2 
alpha2 <- 30
beta2 <- 10

# Generate bimodal data
set.seed(0)
d <- data.frame(y = c(rbeta(100, alpha1, beta1), rbeta(50, alpha2, beta2)))
# correction recommended in cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf, cf first paragraph in section 2, page 3):
n <- length(d$y)
d$yc <- (d$y* (n-1)+0.5)/n

# histogram 
par(mfrow=c(1,2))
hist(d$yc, 50)

# fitting mixtures of beta distributions
uni.modal <- betamix(yc ~ 1 | 1, data = d, k = 1)
bi.modal <- betamix(yc ~ 1 | 1, data = d, k = 2)

# 1) test for bimodality:
lrtest(uni.modal, bi.modal)

# 3) determine for each data point which distribution it belongs
d$group <- (posterior(bi.modal)[,1] <= posterior(bi.modal)[,2])+1
plot(d$y, col=d$group)


# 2) characterize the two underlying distributions
# betamix uses a different parametrization than dbeta (cf cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf). 
# instead of alpha and beta, betareg parametrizes the beta distribution using mu and phi, where
# mu=alpha/(alpha+beta)
# phi= alpha + beta

# ERROR: converting mu and phi to alpha and beta
mu1 <- coef(bi.modal)[1,1]
phi1 <- coef(bi.modal)[1,2]

(a1 <- mu1*phi1) # output: 4.61028
(b1 <- (1-mu1)*phi1) # ouput:-0.7240308


# 4) overaly the density distributions on the histogram

Answer 1

All items are not solved fully correctly in the question. I would recommend the following.

(0) The observations "y" do not need to be corrected as they are between 0 and 1 already. Applying the correction shouldn't create problems but it's not necessary either.

(1) cannot be answered by the likelihood ratio (LR) test. Generally in mixture models, the selection of the number of components cannot be based on the LR test because its regularity assumptions are not fulfilled. Instead, information criteria are often used and "flexmix" upon which betamix() is based offers AIC, BIC, and ICL. So you could choose the best BIC solution among 1, 2, 3 clusters via

library("flexmix")
set.seed(0)
m <- betamix(y ~ 1 | 1, data = d, k = 1:3)

(2) The parameters in betamix() are not mu and phi directly but additionally link functions are employed for both parameters. The defaults are logit and log, respectively. This ensure that the parameters are in their valid ranges (0, 1) and (0, inf), respectively. One could refit the models in both components to get easier access to the links and inverse links etc. However, here it is probably easiest to apply the inverse links by hand:

mu <- plogis(coef(m)[,1])
phi <- exp(coef(m)[,2])

This shows that the means are very different (0.25 and 0.77) while the precisions are rather similar (49.4 and 47.8). Then we can transform back to alpha and beta which gives 12.4, 37.0 and 36.7, 11.1 which is reasonably close to the original parameters in the simulation:

a <- mu * phi
b <- (1 - mu) * phi

(3) The clusters can be extracted using the clusters() function. This simply selects the component with the highest posterior() probability. In this case, the posterior() is really clear-cut, i.e., either close to zero or close to 1.

cl <- clusters(m)

(4) When visualizing the data with histograms, one can either visualize both components separately, i.e., each with its own density function. Or one can draw one joint histogram with the corresponding joint density. The difference is that the latter needs to factor in the different cluster sizes: the prior weights are about 1/3 and 2/3 here. The separate histograms can be drawn like this:

## separate histograms for both clusters
hist(subset(d, cl == 1)$y, breaks = 0:25/25, freq = FALSE,
  col = hcl(0, 50, 80), main = "", xlab = "y", ylim = c(0, 9))

hist(subset(d, cl == 2)$y, breaks = 0:25/25, freq = FALSE,
  col = hcl(240, 50, 80), main = "", xlab = "y", ylim = c(0, 9), add = TRUE)

## lines for fitted densities
ys <- seq(0, 1, by = 0.01)
lines(ys, dbeta(ys, shape1 = a[1], shape2 = b[1]),
  col = hcl(0, 80, 50), lwd = 2)
lines(ys, dbeta(ys, shape1 = a[2], shape2 = b[2]),
  col = hcl(240, 80, 50), lwd = 2)

## lines for corresponding means
abline(v = mu[1], col = hcl(0, 80, 50), lty = 2, lwd = 2)
abline(v = mu[2], col = hcl(240, 80, 50), lty = 2, lwd = 2)

And the joint histogram:

p <- prior(m$flexmix)
hist(d$y, breaks = 0:25/25, freq = FALSE,
  main = "", xlab = "y", ylim = c(0, 4.5))
lines(ys, p[1] * dbeta(ys, shape1 = a[1], shape2 = b[1]) +
  p[2] * dbeta(ys, shape1 = a[2], shape2 = b[2]), lwd = 2)

The resulting figure is included below.