Fit Mixture Of Gamma And Negative Binomial Distributions In R

Question

[FGV]

Hi,

I'm trying to fit a mixture of a gamma and negative binomial distribution to a data set in R. In the end the idea is to get the parameters of the two fitted distributions as well as a likelihood; something like this but, as far as I can tell, this package only does mixtures of the same type of distributions (two normal, two poisson, etc..)

I've looked arround and the R package that looks more promising is flexmix but I cannot get it to work. I've been trying with simpler examples (gaussian and poisson) but no good. So far this is what I've got:

ex1 <- stepFlexmix(y~x,
                   data=data, k=2:5,
                   model=list(FLXMRglm(.~., family="gaussian"),
                              FLXMRglm(.~., family="poisson")),
                   control=list(verb=5, iter=100))

Anyone has any idea on how to get it?

Thanks in advance,

FGV

[Arun]

(may be irrelevant, but out of curiosity) why do you fit with gamma and negative binomial? Negative binomial could already be obtained by modeling $\lambda$ parameter of poisson distribution as a random variable that is gamma distributed. http://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma.E2.80.93Poisson_mixture

[FGV]

I'm analyzing some data that is supposed to be a mixture of two different processes: one has been described as following a negative binomial but the other is a bit more obscure (so I'am trying the gamma). My objective is to try to disentangle both effects.

I know that a neg binomial can be obtained by a poisson and a gamma but am not entirely sure how that can help...

[Zev.Kronenberg]

Is this read depth data?

[FGV]

Yes, the data are NGS read depth

[Arun]

Yes indeed, you're right. It has nothing to do with your question :). I was thinking more in the line of Zev. Reads from RNA-seq data is usually modeled as a NB distribution as well. I thought you had it gamma and NB confused. Sorry about that.

[muhammad.daniyal]

respected sir

I am working on the mixture of burr XII and weibull dstributions. I have to estimate the parameters through R language. can you kindly guide me how I will do it in R. thanx

Answer 1

[Arun]

How about this? you'll need the BB package

# install.packages("BB") # if you already don't have it.
require(BB)

# generate test data to check model efficiency
set.seed(100) 
gamma.nb.mix <- function(n, prob=0.5, shape=2, scale=15, mu=1, sig=0.1) {
    u <- runif(n) 
    out <- apply( as.matrix(u), 1, function(x) ifelse(x<=0.5, rgamma(1, shape=shape, scale=scale), rnorm(1, mu, sig) ) )
    out 
}
out <- gamma.nb.mix(1000)

# maximizing function
f = function(shape, scale, mu, sig, prob) { 
     -sum(log(prob*dgamma(out, shape=shape, scale=scale) + 
       (1-prob)*dnorm(out, mean=mu, sd=sig))) 
}

# arbitrary start (probability parameter is important to remain as close as possible)
start0 <- list("shape"=1.2, "scale"=4, "mu"=0.35, "sig"=0.6, "prob"=0.5) 

fcn <- function(x) f(x[1], x[2], x[3], x[4], x[5])   
res <- spg(par=as.numeric(start0), fn=fcn, lower=c(0, 0, -Inf, 0, 0), upper=c(Inf, Inf, Inf, Inf, 1), control=list( maxit=10000 )  )

After 10k iterations:

Initial model parameters: shape=2, scale=15, mu=1, sig=0.1
Final model parameters (non-converged, but after limit reached):  shape=1.92 scale=15.399  mu=0.999  sigma=0.102  prob=0.5149593

Any good?

[Arun]

In general probability/statistics questions are more appropriate on stats.stackexchange

[Arun]

I just realized that I implemented for gamma and normal distribution. I'll try once again with gamma and NB and report back.