[R] Using betareg package to fit beta mixture with given initial parameters

[Achim Zeileis]

On Wed, 22 Mar 2017, Michael Dayan wrote:

> The method of setting the initial values given lambda, alpha1, etc. should
> not depend on the exact values of lambda, alpha1, etc. in my situation,
> i.e. it does not depend on my data.

Presently, flexmix() that betamix() is built on cannot take the parameters 
directly for initialization. However, it is possible to pass a matrix with 
initial 'cluster' probabilities. This can be easily generated using 
dbeta().

For a concrete example consider the data generated in this discussion on 
SO:

http://stats.stackexchange.com/questions/114959/mixture-of-beta-distributions-full-example

Using that data with random starting values requires 42 iterations until 
convergence:

set.seed(0)
m1 <- betamix(y ~ 1 | 1, data = d, k = 2)
m1

## Call:
## betamix(formula = y ~ 1 | 1, data = d, k = 2)
## 
## Cluster sizes:
##   1   2
##  50 100
## 
## convergence after 42 iterations

Instead we could initialize with the posterior probabilities obtained at 
the observed data (d$y), the true alpha/beta parameters (10; 30 vs. 30; 
10) and the true cluster proportions (2/3 vs. 1/3):

p <- cbind(2/3 * dbeta(d$y, 10, 30), 1/3 * dbeta(d$y, 30, 10))
p <- p/rowSums(p)

This converges after only 2 iterations:

set.seed(0)
m2 <- betamix(y ~ 1 | 1, data = d, k = 2, cluster = p)
m2

## Call:
## betamix(formula = y ~ 1 | 1, data = d, k = 2, cluster = p)
## 
## Cluster sizes:
##   1   2
## 100  50
## 
## convergence after 2 iterations

Up to label switching and small numerical differences, the parameter 
estimates agree. (Of course, these are on the mu/phi scale and not 
alpha/beta as explained in the SO post linked above.)

coef(m1)
##        (Intercept) (phi)_(Intercept)
## Comp.1    1.196286          3.867808
## Comp.2   -1.096487          3.898976
coef(m2)
##        (Intercept) (phi)_(Intercept)
## Comp.1   -1.096487          3.898976
## Comp.2    1.196286          3.867808


> On Mar 22, 2017 04:30, "David Winsemius" <dwinsemius at comcast.net> wrote:
>
>
>> On Mar 21, 2017, at 5:04 AM, Michael Dayan <mdayan.research at gmail.com>
> wrote:
>>
>> Hi,
>>
>> I would like to fit a mixture of two beta distributions with parameters
>> (alpha1, beta1) for the first component, (alpha2, beta2) for the second
>> component, and lambda for the mixing parameter. I also would like to set a
>> maximum of 200 iterations and a tolerance of 1e-08.
>>
>> My question is: how to use the betareg package to run the fit with initial
>> values for the parameters alpha1, beta1, alpha2, beta2 and lambda? I saw
> in
>> the documentation that I would need to use the 'start' option of the
>> betareg function, with start described as "an optional vector with
> starting
>> values for all parameters (including phi)". However I could not find how
> to
>> define this list given my alpha1, beta1, alpha2, beta2 and lambda.
>>
>> The current code I have is:
>> mydata$y <- <my_data>
>> bmix <- betamix(y ~ 1 | 1, data = mydata, k = 2, fstol = 1e-08, maxit =
>> 200, fsmaxit = 200)
>>
>>
>> And I suspect I would need to do something along the lines:
>>
>> initial.vals <- c(?, ?, ?, ?, ?)
>> bmix <- betamix(y ~ 1 | 1, data = mydata, k = 2, fstol = 1e-08, maxit =
>> 200, fsmaxit = 200, control=betareg.control(start=initial.vals)))
>>
>> But I do not know what to use for initial.vals.
>
> If there were sensitivity to data, then wouldn't  that depend on your
> (unprovided) data?
>
>
>>
>> Best wishes,
>>
>> Michael
>>
>>       [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at 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.
>
> David Winsemius
> Alameda, CA, USA
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.
>