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

[Michael Dayan]

Dear Achim,

Thank you for your help, this is exactly what I needed. Your help and
didactic example is very much appreciated.

Best wishes,

Michael

On Wed, Mar 22, 2017 at 11:34 AM, Achim Zeileis <Achim.Zeileis at uibk.ac.at>
wrote:

> 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-b
> eta-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/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>

	[[alternative HTML version deleted]]