4
$\begingroup$

My data is a list of observations and a count for each observation. The data is overdispersed, the mean is ~1,200 and the variance is ~18,000,000. I want to use a negative binomial model to assign p-values for each observation. I tried doing this with a Poisson model (I know the number of trials and the probability of success for each trial) but the p-values became so small for many of the observations that python interpreted the number as 0. Values that were just a couple hundred from the mean were called highly significant because the variance was being underestimated by the model.

Is there an easy to implement way to estimate the parameters for the negative binomial using my data in python? I also need help developing an intuition for what r and p represent in the negative binomial when it is really a mixed Poisson Gamma distribution.

Improve this question
$\endgroup$
6
  • $\begingroup$ After reading this blog post it appears that r=beta and p = alpha/(1-alpha) where alpha and beta are the parameters from the gamma distribution. So I guess I need to estimate alpha and beta somehow. $\endgroup$
    – elsherbini
    Commented Mar 29, 2015 at 21:38
  • $\begingroup$ Not quite sure what you mean by the p-value of a single observation. $\endgroup$
    – Scortchi ♦
    Commented Mar 30, 2015 at 11:07
  • $\begingroup$ Sorry, my data is a list of observations with counts asociated with each one. Something like [100, 42000, 1300, ...]. Once I have a ML estimate of the parameters for the distribution, I want to ask what the chances are of seeing 42000, for instance, by random chance. $\endgroup$
    – elsherbini
    Commented Mar 30, 2015 at 12:26
  • $\begingroup$ Not a p-value then, just a probability under an assumed distribution. And note that small probabilities of individual observations aren't in themselves evidence that the assumed distribution's not a good fit (think of the distribution of winning combinations in a lottery). Using logarithms would let you deal with them more easily. $\endgroup$
    – Scortchi ♦
    Commented Mar 30, 2015 at 12:53
  • $\begingroup$ In python numbers 0.0 < x < 1e323 are interpreted as 0.0, so taking the log yields -infinity. I rejected the poisson values because the mean and variance were three orders of magnitude apart, and knowing the underlying process it is likely inappropriate to identify the counts of so many observations as significant. $\endgroup$
    – elsherbini
    Commented Mar 30, 2015 at 13:01

1 Answer 1

Reset to default
2
$\begingroup$

For maximum-likelihood estimation, you'll need to solve the score equations numerically: see http://en.wikipedia.org/wiki/Negative_binomial_distribution#Maximum_likelihood_estimation. (Or directly maximize the log-likelihood.)

For method-of-moments estimation, following e.g. the parametrization given here, substitute the sample mean & variance for the population mean $\operatorname{E}Y$ & variance $\operatorname{Var}Y$, & solve for the parameters $\mu$ & $\theta$. In this case the estimates are

$$\tilde\mu = \bar y$$

$$ {\tilde\theta}= \frac{\bar y^2}{s_y^2 - \bar y} $$

where $\bar y$ is the sample mean, & $s_y^2$ the sample variance.

Improve this answer
$\endgroup$
6
  • $\begingroup$ This follow-up may betray my ignorance, but why do I want to do method-of-moments estimation after getting the maximum likelihood estimate of the parameters? I started to implement MLE for the parameters but am not having any success. stackoverflow.com/questions/29338152/… $\endgroup$
    – elsherbini
    Commented Mar 30, 2015 at 12:31
  • $\begingroup$ You don't: they're alternative estimators. (I gave the MOM estimator as you asked for something easy.) $\endgroup$
    – Scortchi ♦
    Commented Mar 30, 2015 at 12:46
  • $\begingroup$ Re your MLE implementation: plot the likelihood to see what's going on. In the linked post you seem to be estimating three parameters, which you haven't mentioned here; you also have some rather extreme outliers. $\endgroup$
    – Scortchi ♦
    Commented Mar 30, 2015 at 13:00
  • $\begingroup$ Ah, great idea, I'll plot the likelihood. The third parameter "loc" should be fixed at 0, and I think I just figured out a way to do that. I'm only interested in the negative binomial parameters r and p. $\endgroup$
    – elsherbini
    Commented Mar 30, 2015 at 13:03
  • $\begingroup$ I think one of the main problems I am having is that for the Poisson model, I understood what n and p represented, but I don't have an intuition what r and p mean for the negative binomial when it is being used for these data. $\endgroup$
    – elsherbini
    Commented Mar 30, 2015 at 13:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.