R help - Nov 2016 - Three-component Negative Binomial Mixture: R code

I need a function for R software which computes a mixture of Negative
Binomial distributions with at least three components.

I found on another site the following function "mixnbinom". It works
very
well, but it computes a mixture of only two components:

> mixnbinom=function(y,k1,mu1,k2,mu2,prob,eps=1/100000)
> {
>  new.parms=c(k1,mu1,k2,mu2,prob)
>  err=1
>  iter=1
>  maxiter=100
>  hist(y,probability=T,nclass=30,col="lightgrey",main="The EM
algorithm")
>  xvals=seq(min(y),max(y),1)
>  lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>            (1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="green")
>  while(err>eps){
>      if(iter<=maxiter){
>          lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>                    (1-prob)*dnbinom(xvals,size=k2,mu=mu2),lty=3)
>      }
>      bayes=(prob*dnbinom(y,size=k1,mu=mu1))/((prob*
>      dnbinom(y,size=k1,mu=mu1))+((1-prob)*dnbinom(y,size=k2,mu=mu2)))
>      mu1=sum(bayes*y)/sum(bayes)
>      mu2=sum((1-bayes)*y)/sum((1-bayes))
>      var1=sum(bayes*(y-mu1)^2)/sum(bayes)
>      var2=sum((1-bayes)*(y-mu2)^2)/sum((1-bayes))
>      k1=abs(mu1/((var1/mu1)-1))
>      k2=abs(mu2/((var2/mu2)-1))
>      prob=mean(bayes)
>      old.parms=new.parms
>      new.parms=c(k1,mu1,k2,mu2,prob)
>      err=max(abs((old.parms-new.parms)/new.parms))
>      iter=iter+1
>  }
>  lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>            (1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="red")
>  print(list(k1=k1,mu1=mu1,k2=k2,mu2=mu2,p=prob,iter=iter,err=err))
> }

I would be grateful if someone can modify the previous function to
model a three-component mixture instead of a two-component one.

I also tried to look for a package which does the same job: I have
used the package "mixdist", but I am not able to set up a suitable set
of starting parameters for the function mix (they always converge to
zero). Hereafter, I found the package "DEXUS", but the related
function
does not provide good estimates for the model, even in the event that
I already know what results I have to expect.

Any help is highly appreciated.


Danilo Carit?

-------------------------------------------------------------
Danilo Carit?

PhD Candidate
University of Naples "Parthenope"
Dipartimento di Studi Aziendali e Quantitativi
via G. Parisi, 13, 80132 Napoli - Italy

Opinion only (and therefore ignorable):

Unless you have a lot (?) of data, fitting such a model successfully
is likely to be very challenging. I suggest that you consult a local
statistical expert to see if you might take a different approach.



On Mon, Nov 7, 2016 at 4:42 AM,  <danilo.carita at uniparthenope.it>
wrote:
> I need a function for R software which computes a mixture of Negative
> Binomial distributions with at least three components.
>
> I found on another site the following function "mixnbinom". It
works very
> well, but it computes a mixture of only two components:
>
>
>> mixnbinom=function(y,k1,mu1,k2,mu2,prob,eps=1/100000)
>> {
>>  new.parms=c(k1,mu1,k2,mu2,prob)
>>  err=1
>>  iter=1
>>  maxiter=100
>> 
hist(y,probability=T,nclass=30,col="lightgrey",main="The EM
algorithm")
>>  xvals=seq(min(y),max(y),1)
>>  lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>>           
(1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="green")
>>  while(err>eps){
>>      if(iter<=maxiter){
>>          lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>>                    (1-prob)*dnbinom(xvals,size=k2,mu=mu2),lty=3)
>>      }
>>      bayes=(prob*dnbinom(y,size=k1,mu=mu1))/((prob*
>>      dnbinom(y,size=k1,mu=mu1))+((1-prob)*dnbinom(y,size=k2,mu=mu2)))
>>      mu1=sum(bayes*y)/sum(bayes)
>>      mu2=sum((1-bayes)*y)/sum((1-bayes))
>>      var1=sum(bayes*(y-mu1)^2)/sum(bayes)
>>      var2=sum((1-bayes)*(y-mu2)^2)/sum((1-bayes))
>>      k1=abs(mu1/((var1/mu1)-1))
>>      k2=abs(mu2/((var2/mu2)-1))
>>      prob=mean(bayes)
>>      old.parms=new.parms
>>      new.parms=c(k1,mu1,k2,mu2,prob)
>>      err=max(abs((old.parms-new.parms)/new.parms))
>>      iter=iter+1
>>  }
>>  lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>>            (1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="red")
>>  print(list(k1=k1,mu1=mu1,k2=k2,mu2=mu2,p=prob,iter=iter,err=err))
>> }
>
>
>
> I would be grateful if someone can modify the previous function to
> model a three-component mixture instead of a two-component one.
Isn't that your job?

-- Bert
>
> I also tried to look for a package which does the same job: I have
> used the package "mixdist", but I am not able to set up a
suitable set
> of starting parameters for the function mix (they always converge to
> zero). Hereafter, I found the package "DEXUS", but the related
function
> does not provide good estimates for the model, even in the event that
> I already know what results I have to expect.
>
> Any help is highly appreciated.
>
>
> Danilo Carit?
>
> -------------------------------------------------------------
> Danilo Carit?
>
> PhD Candidate
> University of Naples "Parthenope"
> Dipartimento di Studi Aziendali e Quantitativi
> via G. Parisi, 13, 80132 Napoli - Italy
>
> ______________________________________________
> 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.

On Mon, 7 Nov 2016, danilo.carita at uniparthenope.it wrote:
> I need a function for R software which computes a mixture of Negative
> Binomial distributions with at least three components.
The package "countreg" on R-Forge provides a driver FLXMRnegbin() that
can
be combined with the "flexmix" package (i.e., functions flexmix() and 
stepFlexmix()). The manual page provides some worked illustrations in 
example("FLXMRnegbin", package = "countreg").

Note that the driver is mainly designed for negative binomial _regression_ 
models. But if you just regress on a constant (y ~ 1) you can also get 
negative binomial mixture distributions without covariates.
> I found on another site the following function "mixnbinom". It
works very
> well, but it computes a mixture of only two components:
>
>
>> mixnbinom=function(y,k1,mu1,k2,mu2,prob,eps=1/100000)
>> {
>>  new.parms=c(k1,mu1,k2,mu2,prob)
>>  err=1
>>  iter=1
>>  maxiter=100
>> 
hist(y,probability=T,nclass=30,col="lightgrey",main="The EM
algorithm")
>>  xvals=seq(min(y),max(y),1)
>>  lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>>           
(1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="green")
>>  while(err>eps){
>>      if(iter<=maxiter){
>>          lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>>                    (1-prob)*dnbinom(xvals,size=k2,mu=mu2),lty=3)
>>      }
>>      bayes=(prob*dnbinom(y,size=k1,mu=mu1))/((prob*
>>      dnbinom(y,size=k1,mu=mu1))+((1-prob)*dnbinom(y,size=k2,mu=mu2)))
>>      mu1=sum(bayes*y)/sum(bayes)
>>      mu2=sum((1-bayes)*y)/sum((1-bayes))
>>      var1=sum(bayes*(y-mu1)^2)/sum(bayes)
>>      var2=sum((1-bayes)*(y-mu2)^2)/sum((1-bayes))
>>      k1=abs(mu1/((var1/mu1)-1))
>>      k2=abs(mu2/((var2/mu2)-1))
>>      prob=mean(bayes)
>>      old.parms=new.parms
>>      new.parms=c(k1,mu1,k2,mu2,prob)
>>      err=max(abs((old.parms-new.parms)/new.parms))
>>      iter=iter+1
>>  }
>>  lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+
>>            (1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="red")
>>  print(list(k1=k1,mu1=mu1,k2=k2,mu2=mu2,p=prob,iter=iter,err=err))
>> }
>
>
> I would be grateful if someone can modify the previous function to
> model a three-component mixture instead of a two-component one.
>
> I also tried to look for a package which does the same job: I have
> used the package "mixdist", but I am not able to set up a
suitable set
> of starting parameters for the function mix (they always converge to
> zero). Hereafter, I found the package "DEXUS", but the related
function
> does not provide good estimates for the model, even in the event that
> I already know what results I have to expect.
>
> Any help is highly appreciated.
>
>
> Danilo Carit?
>
> -------------------------------------------------------------
> Danilo Carit?
>
> PhD Candidate
> University of Naples "Parthenope"
> Dipartimento di Studi Aziendali e Quantitativi
> via G. Parisi, 13, 80132 Napoli - Italy
>
> ______________________________________________
> 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.

I tried the function flexmix() with the driver FLXMRnegbin() with two  
components first, in order to compare its results with those provided  
by my function mixnbinom(). In particular, I ran the following code:

> fm0 <- flexmix(y ~ 1, data = data.frame(y), k = 2, model =
FLXMRnegbin())

where "y" is my vector of counts. The previous function provided me  
the following parameters:

>                    Comp.1   Comp.2
> coef.(Intercept) 1.2746536 1.788578
> theta            0.1418201 5.028766

with priors 0.342874 and 0.657126, respectively. I assume that the  
coefficients "Intercept" represent the two means of the model (mu1 and
mu2), while the "theta" coefficients are the size parameters (size1  
and size2).

Unfortunately, unlike my function mixnbinom(), the model computed with  
flexmix() did not provide a good fit to my data (p-value ~0).

Is there something wrong in the process above?



Achim Zeileis <Achim.Zeileis at uibk.ac.at> ha scritto:
> On Mon, 7 Nov 2016, danilo.carita at uniparthenope.it wrote:
>
>> I need a function for R software which computes a mixture of Negative
>> Binomial distributions with at least three components.
>
> The package "countreg" on R-Forge provides a driver FLXMRnegbin()
> that can be combined with the "flexmix" package (i.e., functions
> flexmix() and stepFlexmix()). The manual page provides some worked  
> illustrations in example("FLXMRnegbin", package =
"countreg").
>
> Note that the driver is mainly designed for negative binomial  
> _regression_ models. But if you just regress on a constant (y ~ 1)  
> you can also get negative binomial mixture distributions without  
> covariates.

-------------------------------------------------------------
Danilo Carit?

PhD Candidate
University of Naples "Parthenope"
Dipartimento di Studi Aziendali e Quantitativi
via G. Parisi, 13, 80132 Napoli - Italy