Counting samples seem to not be Poisson distributed, need sanity check

Question

I have an exercise where I have to use Poisson one-way classification / Regression of some data. The data I have is a set of 120 samples grouped the following labels A, B, C, D, E, and F. For each group, there are 20 samples (or 20 repetitions) with a count value. Now that is all good and from what I can tell is well suited for the assumption that it may fit a Poisson distribution.

However, as I understand it one of the properties of a random variable $Y$ that follows $$Y \sim Po(\lambda)$$ Then it follows that $$E(Y) = V(Y)$$ But when I calculate the means (expected) and the variance for the data according to the grouping: Then I get

|          |         A |         B |         C |         D |          E |         F |
|----------+-----------+-----------+-----------+-----------+------------+-----------|
|----------+-----------+-----------+-----------+-----------+------------+-----------|
| Mean     |      4.90 |      9.45 |      8.65 |      1.45 |      18.35 |      0.80 |
| Variance | 9.8842105 | 6.4710526 | 7.3973684 | 1.3131579 | 15.6078947 | 0.5894737 |
|----------+-----------+-----------+-----------+-----------+------------+-----------|

So from what I can tell $$E(Y) \neq V(Y)$$ and just to show how I computed this using R:

( Means <- tapply(D$NumberPGrains, D$Era, mean) )
( Variances <- tapply(D$NumberPGrains, D$Era, var) )

This means, from my understanding that the data is not Poisson distributed. So my question is: Am I wrong can this still be Poisson distributed, is there something I have missed?

And just to clarify, the exercise literally states to follow Poisson one-way classification (the title of the exercise: "Question 3 -Poisson one-way classification model"), but right now I have hard time seeing the purpose of that.

Answer 1

A true Poisson distribution will have its mean exactly equal to its variance. For a sampling of a Poisson distribution, however, there will be some deviation - with only 20 samples, it's unlikely that you'd see the mean and variance of the sample be exactly equal. For the most part, you seem to have a strong correlation between the mean and variance, which is good. You could also find the confidence intervals around your parameter estimates, to take a hypothesis testing approach to determine whether your mean and variance estimates are really statistically different from one another. With a very large sample size, you'll have very good estimates which should be very nearly equal if the data is indeed Poisson distributed, but for lower sample sizes, your estimates won't be as good, so some numerical differences between the mean and variance are expected.

Answer 2

+1 to Nuclear Wang's answer.

It's always good to simulate a little to get a feeling for the randomness involved in a situation like this. Here is a little something I hacked together, which may be enlightening.

I'll take the means you observed as the true means of your underlying data generating process. Then I'll assume that the data are actually Poisson distributed, and generate simulations using these assumptions. Then I'll calculate the variance for such simulations. I'll do this, say, fifty times, so I end up with fifty different variances we might observe under the null hypothesis of Poisson distributions with the means you observed. Finally, I'll plot these variances (as dot coulds, jittered horizontally), and also the variances you actually observed (as horizontal red lines):

In panels B-F, the observed variances are pretty much in the middle of the point clouds of simulated variances, so they seem to be quite consistent with a Poisson assumption. In panel A, the variance you observe is way higher than all the variances we observed, so this one does not appear consistent with a Poisson assumption. (We could call this an informal hypothesis test and say that we got $p=\frac{1}{50}=0.02$, though this doesn't account for the uncertainty in the mean.)

R code:

n_sims <- 50
n_sample <- 20
means <- c(4.90, 9.45, 8.65, 1.45, 18.35, 0.80)
vars_obs <- c(9.8842105, 6.4710526, 7.3973684, 1.3131579, 15.6078947, 0.5894737)

vars <- matrix(NA,nrow=n_sims,ncol=length(means),
    dimnames=list(NULL,LETTERS[seq_along(means)]))
for ( ii in 1:n_sims ) {
    set.seed(ii)    # for reproducibility
    vars[ii,] <- sapply(means,function(mm)var(rpois(n_sample,mm)))
}

opar <- par(mfrow=c(2,3))
    for ( jj in seq_along(means) ) {
        plot(rnorm(n_sims),vars[,jj],pch=19,main=colnames(vars)[jj],
            ylim=range(c(vars[,jj],vars_obs[jj])),xaxt="n",las=1,ylab="",xlab="")
        abline(h=vars_obs[jj],col="red",lwd=2)
    }
par(opar)

Answer 3

One can formally test the fit of the data to the theoretical poisson distribution using a chi-squared goodness of fit test, via the goodfit() function included in the vcd package.

• Null hypothesis significance test: If p < .05 or whatever
  threshold we might choose in advance, we will reject the hypothesis
  that x comes from a Poisson process.

• More flexible approach: we simply assess the degree of evidence
  against the null based on the size of the p-value, where smaller
  p = stronger evidence.

  require(vcd)
  
  x=0:5
  
  freqx=c(7,12,7,3,1,1) # Counts at each level of x
  gf = goodfit (x, type= "poisson", method= "ML")
  

summary(gf)

    Goodness-of-fit test for poisson distribution

                      X^2 Likelihood Ratio = 1.91, df = 4, P(> X^2) = 0.75


plot (gf, main= "x vs Poisson distribution")

In plots such as these from vcd, bars lifted above the x-axis reflect negative residuals (too few occurrences to fit the theoretical distribution); bars that extend below, positive residuals (too many). This may seem less than intuitive, and the plot can't be modified in some of the usual ways (text size relative to bar size, etc), and so the following approach may be more helpful.

require(fitdistrplus)
poisson = fitdist(x, 'pois', method = 'mle')
print(poisson)  # Identify the best value of lambda (which indicates both mean and variance) to insert below
dist = dpois(0:5, lambda = 2.5) * sum(freqx) # assign lambda from the fit above
df = as.data.frame(dist) 
windows()
# Then we plot first the observed distribution, as vertical lines, and then the theoretical distribution as the curve
plot(x, freqx,  type='h', lwd=2, main = 'Curve = Fitted Poisson Values',
  cex.axis=1.2, cex.lab=1.2)
lines(x, df$dist, col = 'red', lwd=3)