

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 1e08.
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 = 1e08, 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 = 1e08, maxit =
200, fsmaxit = 200, control=betareg.control(start=initial.vals)))
But I do not know what to use for initial.vals.
Best wishes,
Michael
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


> On Mar 21, 2017, at 5:04 AM, Michael Dayan < [hidden email]> 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 1e08.
>
> 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 = 1e08, 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 = 1e08, 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?
David Winsemius
Alameda, CA, USA
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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.
On Mar 22, 2017 04:30, "David Winsemius" < [hidden email]> wrote:
> On Mar 21, 2017, at 5:04 AM, Michael Dayan < [hidden email]>
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 1e08.
>
> 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 = 1e08, 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 = 1e08, 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?
postingguide.html
> and provide commented, minimal, selfcontained, reproducible code.
David Winsemius
Alameda, CA, USA
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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/mixtureofbetadistributionsfullexampleUsing 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" < [hidden email]> wrote:
>
>
>> On Mar 21, 2017, at 5:04 AM, Michael Dayan < [hidden email]>
> 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 1e08.
>>
>> 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 = 1e08, 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 = 1e08, 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]]
>>
>> ______________________________________________
>> [hidden email] mailing list  To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide http://www.Rproject.org/> postingguide.html
>> and provide commented, minimal, selfcontained, reproducible code.
>
> David Winsemius
> Alameda, CA, USA
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list  To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
>
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Tree dbh haut Btot
1 35.00 18.90 0.357
2 25.00 16.60 0.214
3 23.00 19.50 0.173
4 13.50 15.60 0.060
5 20.00 18.80 0.134
6 23.00 17.40 0.137
7 29.00 19.90 0.428
8 17.60 18.20 0.100
9 31.00 25.30 0.514
10 26.00 23.50 0.273
11 13.00 13.00 0.031
12 32.00 20.70 0.356
13 28.00 28.50 0.349
14 15.00 18.20 0.068
15 19.00 14.90 0.110
16 33.00 14.90 0.260
17 24.00 19.10 0.160
18 22.00 19.20 0.204
19 39.00 25.20 0.724
20 30.00 26.60 0.386
Hello dear all:
I am using above data to run the following code;
library(nlme)
start < coef(lm(Btot~I(dbh**2*haut),data=dat))
names(start) < c("a","b")
model1 <(nlme (Btot~a+b*dbh**2*haut, data=cbind(dat, g="a"), fixed=a+b~1,
start=start,
groups=~g, weights=varPower(form=~dbh)))
I get regression parameters with the intercept being nonsignificant.
Therefore, I run the following code to obtain an equation without
intercept..
summary(nlme(Btot~b*dbh**2*haut, data=cbind(dat,g="a"), fixed=b~1,
start=start["b"], groups=~g, weights=varPower(form=~dbh)))
When I do run the last code, I get the following error:
Error in nlme.formula(Btot ~ b * dbh**2 * haut, data = cbind(dat, g = "a"),
:
step halving factor reduced below minimum in PNLS step
Can someone help??
Best regards,
Santiago
On Wed, Mar 22, 2017 at 2:26 AM, Michael Dayan < [hidden email]>
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.
>
> On Mar 22, 2017 04:30, "David Winsemius" < [hidden email]> wrote:
>
>
> > On Mar 21, 2017, at 5:04 AM, Michael Dayan < [hidden email]>
> 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 1e08.
> >
> > 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 = 1e08, 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 = 1e08, 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]]
> >
> > ______________________________________________
> > [hidden email] mailing list  To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/rhelp> > PLEASE do read the posting guide http://www.Rproject.org/> postingguide.html
> > and provide commented, minimal, selfcontained, reproducible code.
>
> David Winsemius
> Alameda, CA, USA
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list  To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/> postingguide.html
> and provide commented, minimal, selfcontained, reproducible code.
>
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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 < [hidden email]>
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/mixtureofb> etadistributionsfullexample
>
> 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" < [hidden email]> wrote:
>>
>>
>> On Mar 21, 2017, at 5:04 AM, Michael Dayan < [hidden email]>
>>>
>> 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 1e08.
>>>
>>> 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 = 1e08, 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 = 1e08, 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]]
>>>
>>> ______________________________________________
>>> [hidden email] mailing list  To UNSUBSCRIBE and more, see
>>> https://stat.ethz.ch/mailman/listinfo/rhelp>>> PLEASE do read the posting guide http://www.Rproject.org/>>>
>> postingguide.html
>>
>>> and provide commented, minimal, selfcontained, reproducible code.
>>>
>>
>> David Winsemius
>> Alameda, CA, USA
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> [hidden email] mailing list  To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide http://www.Rproject.org/posti>> ngguide.html
>> and provide commented, minimal, selfcontained, reproducible code.
>>
>>
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

