6
$\begingroup$

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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Obviously you have accomplished (3): the points are labeled by color in the right hand plot. What, then, is your question? What goes wrong? $\endgroup$ – whuber Sep 10 '14 at 15:29
  • $\begingroup$ How does one compute shape1 (alpha) and shape2 (beta) parameters from the betamix parameters (mu, phi)? The way I did it leads to a negative beta. How to overlay the fitted density distributions on the histogram? $\endgroup$ – user3464896 Sep 11 '14 at 6:03
11
$\begingroup$

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.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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