# Fitting a probability distribution to zero inflated data in R

I am trying to learn how to fit a probability distribution to a vector of data, using the program R, but there are a lot of potential probability distributions to use! So my question is, how do I find the best distribution for my data, and how do I prove that I have picked the right distribution? Can I acquire AIC values for a whole set of different distributions?

The data are observational count data of bees visiting flowers. Each species has a certain number of visits, hence the differing frequencies. The goal is to find the best distribution to describe the bee visitation, show that I have selected the right one, and then use that distribution to sample from randomly for a set of simulations.

Here is what the data looks like, it is a vector of count observations. It is zero inflated, with a long tailed distribution (maybe zero-inflated negative binomial?).

i.vec=c(0,63,1,4,1,44,2,2,1,0,1,0,0,0,0,1,0,0,3,0,0,2,0,0,0,0,0,2,0,0,0,0,
0,0,0,0,0,0,0,0,6,1,11,1,1,0,0,0,2)


And here are some basic parameters that I have calculated. I am using standard deviation for sigma, and phi is the proportion of zeroes in the data.

m=mean(i.vec)
#[1] 3.040816
sig=sd(i.vec)
#[1] 10.86078
tab<-table(i.vec)
zero.prop<-as.numeric(tab[1])/sum(as.numeric(tab))
#[1] 0.6122449


As you can see, the standard deviation is much greater than the mean, and I have a very high proportion of zeroes.

-
It would be helpful to have some context: where do the data come from? Why do you want to fit a distribution? To me, it's the 11, 44, and 63, rather than the 0's, that are unusual. –  Karl Aug 26 '11 at 1:59
I will add in some information to my question. –  Laura Aug 28 '11 at 22:45

You can use Vuong test in pscl package to compare non-nested models. Here is an example

> m1 <- zeroinfl(i.vec ~ 1 | 1, dist = "negbin")
> summary(m1)

Call:
zeroinfl(formula = i.vec ~ 1 | 1, dist = "negbin")

Pearson residuals:
Min      1Q  Median      3Q     Max
-0.3730 -0.3730 -0.3730 -0.2503  7.3544

Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)   1.1122     0.3831   2.903  0.00369 **
Log(theta)   -1.9256     0.2839  -6.784 1.17e-11 ***

Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -9.815     96.462  -0.102    0.919
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Theta = 0.1458
Number of iterations in BFGS optimization: 579
Log-likelihood: -80.51 on 3 Df

> m2 <- zeroinfl(i.vec ~ 1 | 1, dist = "poisson")
> summary(m2)

Call:
zeroinfl(formula = i.vec ~ 1 | 1, dist = "poisson")

Pearson residuals:
Min      1Q  Median      3Q     Max
-0.7242 -0.7242 -0.7242 -0.4860 14.2795

Count model coefficients (poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.05911    0.08205    25.1   <2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.4561     0.2933   1.555     0.12
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 11
Log-likelihood: -233.7 on 2 Df

> vuong(m1, m2)
Vuong Non-Nested Hypothesis Test-Statistic: 1.946095
(test-statistic is asymptotically distributed N(0,1) under the
null that the models are indistinguishible)
in this case:
model1 > model2, with p-value 0.02582165


Vuong test also suggests that the zero-inflated negative binomial provides a better fit to your data compared to the ordinary negative binomial (not shown here, but you can fit both models and compare them).

-
Thank you, this is so helpful! –  Laura Aug 29 '11 at 6:53

I don't think you necessarily need to inflate zeros.

Your data seem quite consistent with a negative binomial:

> library(MASS)
> table(rnegbin(49,mu=3.1,theta=0.075))

0  1  2  3  5 18 20 21 31 61
36  4  2  1  1  1  1  1  1  1
> table(i.vec)
i.vec
0  1  2  3  4  6 11 44 63
30  8  5  1  1  1  1  1  1


So lets get some estimates:

> mean(i.vec)
[1] 3.040816
> theta.ml(i.vec,3.041)
[1] 0.145777
attr(,"SE")
[1] 0.04136887


So let's just take a look at mu = 3.041 and theta at say 0.14 (there's lots of uncertainty in theta here):

Here's three random samples from that distribution:

> table(rnegbin(49,mu=3.041,theta=0.14))

0  1  2  3  4  5  6  7  8 18 49
30  5  1  3  2  2  1  2  1  1  1
> table(rnegbin(49,mu=3.041,theta=0.14))

0  1  2  3  4  7  9 15 29 31 47 56
33  2  4  1  1  1  1  2  1  1  1  1
> table(rnegbin(49,mu=3.041,theta=0.14))

0  1  2  3  4  7  9 12 16 48 66
36  4  1  1  1  1  1  1  1  1  1


They look similar enough to your data. Negative binomial seems at least plausible.

You may like to play with the functions in MASS

-
Thank you! It won't let me vote up, but I appreciate your help. –  Laura Aug 28 '11 at 22:47