# Simple comparison of two Poisson means in R?

I have two samples and I would like to determine whether the difference between them is statistically significant or not:

Because this is clearly data of small counts that cannot be approximated by a normal distribution, I don't think a t-test is appropriate. Instead, I believe I may assume that these counts follow (two distinct) Poisson distributions. As stated here, each of these Poisson distributions is best summarized by the respective sample means, as follows:

$$\overline{y} = \frac{1}{N} \sum_{i = 1}^n{y_i}$$

This obviously yields non-integer means (0.4 and 2.6, respectively), so I cannot use functions like poisson.test() or any of the functions from the exactci library directly.

I am aware of this answer regarding the C-test and E-test - but is there a straightforward implementation in R that would do this?

The only reasonable way to do this using poisson.test() I can think of would be to sum the rates for each condition and use the numbers of counts as the respective T parameters:

# Count data for each respective Condition
Cond1 <- c(0, 0, 0, 1, 1)
Cond2 <- c(1, 2, 3, 3, 4)

poisson.test(
x = c(sum(Cond1), sum(Cond2)),
T = c(length(Cond1), length(Cond2)),
alternative = "two.sided"
)


which yields the following result:

    Comparison of Poisson rates

data:  c(sum(Cond1), sum(Cond2)) time base: c(length(Cond1), length(Cond2))
count1 = 2, expected count1 = 7.5, p-value = 0.007385
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
0.01685531 0.67955077
sample estimates:
rate ratio
0.1538462


Is this valid or is there a better way to do this in R?

• Do you know how to do a Poisson regression?
– Dave
Oct 7, 2021 at 12:50
• It's valid: see the help page for poisson.test for specifics.
– whuber
Oct 7, 2021 at 13:23
• @whuber is right. The poisson.test function is valid here and gives you the p-value for the C-test. In general, this method is by far more powerful than the GLM approach below for testing the difference between two Poisson rates with small sample sizes like this. Oct 8, 2021 at 7:42

This is what I would do which I believe it is what @Dave is hinting at in comment.

Fit GLM with Poisson familiy distribution and see the effect of "condition" on the mean of each group.

Prepare the dataset:

Cond1 <- c(0, 0, 0, 1, 1)
Cond2 <- c(1, 2, 3, 3, 4)

dat <- data.frame(
cond= c(rep('C1', length(Cond1)), rep('C2', length(Cond2))),
count= c(Cond1, Cond2)
)


Fit the model and assess the significance of "condition":

fit <- glm(count ~ cond, data= dat, family= 'poisson')
summary(fit)

Call:
glm(formula = count ~ cond, family = "poisson", data = dat)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-1.13533  -0.89443  -0.07296   0.65703   0.80391

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.9163     0.7071  -1.296   0.1950
condC2        1.8718     0.7595   2.464   0.0137 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 14.8823  on 9  degrees of freedom
Residual deviance:  5.8682  on 8  degrees of freedom
AIC: 27.731

Number of Fisher Scoring iterations: 5


So the mean of Cond 1 is estimated as exp(-0.9163) = 0.4, the mean of Cond 2 is exp(-0.9163 + 1.8718) = 2.6 and the difference has p = 0.0137 for the null hypothesis of being zero.

• Thank you, this makes sense and is super useful - especially when comparing more than 2 conditions. I'm just wondering then why poisson.test() gives a different p-value? Is poisson.test() simply wrong to use for this scenario, or does it estimate the p-value somehow differently? Oct 7, 2021 at 13:29
• @Brunox13 I wasn't even aware of poisson.test until now so I can't answer about the difference (and I'm curious myslef). However, since poisson.test takes in input only the sums and the sizes of the vectors, I'm tempted to say that the GLM approach is preferable because it makes better use of the information available. Oct 7, 2021 at 13:34
• @Sextus the glm is a log-link by default, so the glm tests $\log(\lambda_2)-\log(\lambda_1)=0$, which is the same thing as poisson.test. The glm is of course relying on an asymptotic approximation. Oct 8, 2021 at 0:07