Statistical testing: do count data come from the same distribution?

Question

The data I am dealing with are groups of counts, $n_i, i=1..K$. More than a half of these counts are zeros. The null hypothesis is that all the counts come from the same distribution, e.g., Poisson with parameter $\lambda$ $$P(n|\lambda)=\frac{\lambda^n}{n!}e^{-\lambda},$$ in which case I can perform the estimate of the parameter as the mean over all the counts, $$\hat{\lambda}=\frac{1}{K}\sum_{i=1}^Kn_i. $$

There might be however situations where the zero and non-zero counts are generated by different distributions (possibly more than two). I need a statistical test to identify such cases.

Clarification: the problem is not to test whether the distribution is Poissonian, but whether all counts come from the distribution or not. Zeros may be due to $\lambda$ being small... or because they are permanently zero. (I realized after the discussion in the comments, that the initial formulation of my questions is ambiguous.)

Answer 1

You can test the null hypothesis that the data follows a Poisson distribution against the alternative of a zero-inflated Poisson using for example the glmmTMB R-package. The test relies on the likelihood ratio test statistic being asymptotically chi-square with one degree of freedom.

However, unless you also allow zero-deflation under your alternative hypothesis, the null hypothesis of no zero-inflation is on the boundary of the parameter space, and the likelihood ratio statistic is then asymptotically distributed as an equal-weights mixture of chi-squares with zero and one degree of freedom (Stram and Lee 1994) so a better approximate $p$-value would be half of what is computed below.

# Simulated data from a zero-inflated poisson
set.seed(1)
y <- rpois(100, lambda = 3)
y[1:10] <- 0
data <- data.frame(y)

# Testing Poisson against zero-inflated Poisson relying on asymptotic distribution of likelihood ratio statistic
library(glmmTMB)
mod0 <- glmmTMB(y ~ 1, family=poisson, data)
mod1 <- update(mod0, ziformula = ~ 1)
anova(mod0, mod1)
#> Data: data
#> Models:
#> mod0: y ~ 1, zi=~0, disp=~1
#> mod1: y ~ 1, zi=~1, disp=~1
#>      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
#> mod0  1 383.93 386.53 -190.96   381.93                           
#> mod1  2 380.88 386.09 -188.44   376.88 5.0513      1    0.02461 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1