R help - Mar 2017 - Using betareg package to fit beta mixture with given initial parameters

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.

Best wishes,

Michael

	[[alternative HTML version deleted]]

> 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

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" <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]]

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.
>

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 non-significant.
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 <mdayan.research at
gmail.com>
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" <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.
>
	[[alternative HTML version deleted]]