I have a vector of measured proportions and want to
- test if these proportions follow a bimodal distribution
- characterize the two underlying distributions
- determine for each data point which distribution it belongs to
- 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