[R] Mixed Beta Disrubutions

[Bob Aronoff]

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")
}

____