Code described in blog post at http://datavoreconsulting.com/programming-tips/count-data-glms-choosing-poisson-negative-binomial-zero-inflated-poisson/

CountGLM.R

n = 100                             # No. of samples per treatment
mean.A = 10                         # Mean count for treatment A
mean.B = 5                          # Mean count for treatment B
nd = data.frame(Trt = c("A", "B"))  # Used in predict( ) function

# Generate example data based on Poisson distribution (mean = variance)
data.pois = data.frame(Trt = c(rep("A", n), rep("B", n)), 
                       Response = c(rpois(n, mean.A), rpois(n, mean.B))
                       )

# Generalized linear model with Poisson distribution of errors
model.pois = glm(Response ~ Trt, data = data.pois, family = poisson)
summary(model.pois)
# Test for goodness of fit of model; compare residual deviance with df (from summary output)
# If test is significant, model fit is not good
1 - pchisq(summary(model.pois)$deviance, 
           summary(model.pois)$df.residual
           )

# Predict mean counts and standard errors
cbind(nd, 
      Mean = predict(model.pois, newdata=nd, type="response"), 
      SE = predict(model.pois, newdata=nd, type="response", se.fit=T)$se.fit
      )

#--------------------------------------------------------------------------
# Generate example data based on negative binomial distribution (mean < variance)
library(MASS)
data.nb = data.frame(Trt = c(rep("A", n), rep("B", n)), 
                     Response=c(rnegbin(n, mean.A, 5), rnegbin(n, mean.B, 5))
                     )

# Generalized linear model with Poisson distribution of errors
model.pois.2 = glm(Response ~ Trt, data = data.nb, family = poisson)
summary(model.pois.2)
# Test for goodness of fit of model; compare residual deviance with df (from summary output)
1 - pchisq(summary(model.pois.2)$deviance,
           summary(model.pois.2)$df.residual
           )

# Predict mean counts and standard errors
cbind(nd, 
      Mean = predict(model.pois.2, newdata = nd, type = "response"), 
      SE = predict(model.pois.2, newdata = nd, type = "response", se.fit = T)$se.fit
      )

# Generalized linear model with negative binomial distribution of errors
model.nb = glm.nb(Response ~ Trt, data = data.nb)
summary(model.nb)
# Test for goodness of fit of model; compare residual deviance with df (from summary output)
1 - pchisq(summary(model.nb)$deviance,
           summary(model.nb)$df.residual
           )

# Predict mean counts and standard errors
cbind(nd, 
      Mean = predict(model.nb, newdata = nd, type = "response"), 
      SE = predict(model.nb, newdata = nd, type="response", se.fit = T)$se.fit
      )

#--------------------------------------------------------------------------
# Generate example data based on zero-inflated Poisson distribution
library(VGAM)
data.zip = data.frame(Trt = c(rep("A", n), rep("B", n)), 
                      Response = c(rzipois(n, mean.A, 0.2), rzipois(n, mean.B, 0.2))
                      )

# Generalized linear model with Poisson distribution of errors
model.pois.3 = glm(Response ~ Trt, data = data.zip, family = poisson)
summary(model.pois.3)
# Test for goodness of fit of model; compare residual deviance with df (from summary output)
1 - pchisq(summary(model.pois.3)$deviance,
           summary(model.pois.3)$df.residual
           )

# Predict mean counts and standard errors
cbind(nd, 
      Mean = predict(model.pois.3, newdata = nd, type = "response"), 
      SE = predict(model.pois.3, newdata = nd, type = "response", se.fit = T)$se.fit
      )

# Generalized linear model with negative binomial distribution of errors
model.nb.2 = glm.nb(Response ~ Trt, data = data.zip)
summary(model.nb.2)
# Test for goodness of fit of model; compare residual deviance with df (from summary output)
1 - pchisq(summary(model.nb.2)$deviance,
           summary(model.nb.2)$df.residual
           )

# Predict mean counts and standard errors
cbind(nd, 
      Mean = predict(model.nb.2, newdata = nd, type = "response"), 
      SE = predict(model.nb.2, newdata = nd, type = "response", se.fit = T)$se.fit
      )

# Zero-inflated Poisson model
library(pscl)
# Test for treatment effect on counts
model.zip = zeroinfl(Response ~ Trt|1, data = data.zip)
summary(model.zip)

# Predict mean counts and zeros
cbind(nd, 
      Count = predict(model.zip, newdata = nd, type = "count"),
      Zero = predict(model.zip, newdata = nd, type = "zero")
      )

# Test for treatment effect on counts and zeros
model.zip.2 = zeroinfl(Response ~ Trt|Trt, data = data.zip)
summary(model.zip.2)

# Predict mean counts and zeros
cbind(nd, 
      Count = predict(model.zip.2, newdata = nd, type = "count"),
      Zero = predict(model.zip.2, newdata = nd, type = "zero")
      )

# Zero-inflated negative binomial model
model.zip.3 = zeroinfl(Response ~ Trt|1, data = data.zip, dist = "negbin")
summary(model.zip.3)