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