14

I have a vector of count data that is strongly over dispersed and zero inflated.

The vector looks like this:

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)
m=mean(i.vec)
# 3.040816
sig=sd(i.vec)
# 10.86078

I would like to fit a distribution to this, which I strongly suspect will be a zero inflated poisson (ZIP). But I need to perform a significance test to demonstrate that a ZIP distribution fits the data.

If I had a normal distribution, I could do a chi square goodness of fit test using the function goodfit() in the package vcd, but I don't know of any tests that I can perform for zero inflated data.

Tetsujin no Oni
  • 7,300
  • 2
  • 29
  • 46
Laura
  • 679
  • 2
  • 5
  • 14

1 Answers1

22

Here is one approach

# LOAD LIBRARIES
library(fitdistrplus)    # fits distributions using maximum likelihood
library(gamlss)          # defines pdf, cdf of ZIP


# FIT DISTRIBUTION (mu = mean of poisson, sigma = P(X = 0)
fit_zip = fitdist(i.vec, 'ZIP', start = list(mu = 2, sigma = 0.5))

# VISUALIZE TEST AND COMPUTE GOODNESS OF FIT    
plot(fit_zip)
gofstat(fit_zip, print.test = T)

Based on this, it does not look like ZIP is a good fit.

Ramnath
  • 54,439
  • 16
  • 125
  • 152
  • 7
    +1 ; But it does look reasonably acceptable with a negative binomial fit: fit_zip2 <- fitdist(i.vec, 'nbinom', start = list(mu = 3, size = 0.1)) – IRTFM Aug 23 '11 at 14:43
  • 4
    it does look good. i would also try the zero-inflated negative binomial distribution which is available as `zinegbin` in the `VGAM` package. – Ramnath Aug 23 '11 at 15:15
  • 4
    Thank you so much! How did you select your mu and sigma? – Laura Aug 23 '11 at 22:44
  • Also, did you mean to write "print.test=T"? – Laura Aug 23 '11 at 22:49
  • 4
    they are just starting values for the optimization algorithm which i filled arbitrarily. instead you can use `mu = mean(i.vec[i.vec > 0])` and `p` equal to the proportion of zeros in your sample as a starting point. and yes i mean to write `print.test = T` – Ramnath Aug 24 '11 at 01:12