3
$\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.

$\endgroup$
  • $\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 Mar 29 '15 at 21:38
  • $\begingroup$ Not quite sure what you mean by the p-value of a single observation. $\endgroup$ – Scortchi - Reinstate Monica Mar 30 '15 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 Mar 30 '15 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 - Reinstate Monica Mar 30 '15 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 Mar 30 '15 at 13:01
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.

| cite | improve this answer | |
$\endgroup$
  • $\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 Mar 30 '15 at 12:31
  • $\begingroup$ You don't: they're alternative estimators. (I gave the MOM estimator as you asked for something easy.) $\endgroup$ – Scortchi - Reinstate Monica Mar 30 '15 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 - Reinstate Monica Mar 30 '15 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 Mar 30 '15 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 Mar 30 '15 at 13:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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