I am working to understand the same issues with my datasets.
Adapting Dr Zeileis' posting I have written a humble function.
It creates a two cluster beta mix on a vector of data.
It seems to be working well on my datasets.
You are welcome to try on yours.
regards,
Bob
_____
bi.modal.beta<-function (vals){
# vals is a vector data points in range (0,1). There should not be any 0,1,
or NA in vals
library(betareg)
library(flexmix)
d <- data.frame(y = vals)
#create bimodal beta model
m <- betamix(y ~ 1 | 1, data = d, k = 2)
#extract beta parameters (mean and precision) with anti-link functions
mu <- plogis(coef(m)[,1])
phi <- exp(coef(m)[,2])
#convert mean & precision to alpha & beta
a <- mu * phi
b <- (1 - mu) * phi
#report parameters
print("alpha:")
print( a)
print("beta:")
print(b)
#plot in order to inspect result
ys <- seq(0, 1, by = 0.01)
p <- prior(m$flexmix)
# p is equivalent to the percentage of each distribution
hist(d$y, breaks = 0:25/25, freq = FALSE,
main = "Bimodal Betamix", xlab = "values")
lines(ys, p[1] * dbeta(ys, shape1 = a[1], shape2 = b[1]) , lwd = 2, col="red")
lines(ys, p[2] * dbeta(ys, shape1 = a[2], shape2 = b[2]) , lwd = 2,col="blue")
}
____
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
