Unrealistic coefficients when comparing two count samples with glm poisson

Question

I would like to test the significance of the difference in mean between two independent count samples. I'm doing this with a GLM poisson in R, as shown in the code below:

a=c(0,0,0,0,0,0,0,0,0,0)
b=c(1,2,0,1,1,2,0,1,0,2)

c=data.frame(sp=c(a,b),grp=c(rep('A',10),rep('B',10)))

summary(glm(sp~grp,data=c,family=poisson))

Call:
glm(formula = sp ~ grp, family = poisson, data = c)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.41421  -0.00006  -0.00006   0.00000   0.87897  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    -20.3     4914.8  -0.004    0.997
grpB            20.3     4914.8   0.004    0.997

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 22.1807  on 19  degrees of freedom
Residual deviance:  8.3178  on 18  degrees of freedom
AIC: 28.159

Number of Fisher Scoring iterations: 18

As you can see, the coefficient values are not reflecting reality. I noticed that this happened because group "A" has mean=0. In this way, I would like to know if there is any way to fix this problem in glm, or if there is any other better method to test my hypothesis.

Answer 1

The expected value predicted by Poisson regression is

$$ E[y|X,\beta] = e^{\beta_0 + \beta_1 X} $$

assuming we're using the standard log link function. So your model predicts:

> beta <- c(-20.3, 20.3)
> exp(beta[1] + beta[2] * 0)
[1] 1.52694e-09
> exp(beta[1] + beta[2] * 1)
[1] 1
> mean(a)
[1] 0
> mean(b)
[1] 1

I would say that there's nothing wrong with the predictions. The only issue is that you are assuming Poisson distribution for a group that has zero mean, while the distribution is parametrized by $\lambda \ge 0$.

Answer 2

If willing to do an alternative test, here is code to do an E-test for the difference in Poisson means:

# Forgive me, will make putting this in a package on my ToDo list
source('https://dl.dropboxusercontent.com/s/r4zmlxmi41tqsts/ETest.R?dl=0')
poisson.etest(0,10,10,10)

Which returns a p-value of 0.0005782135.

If you want/need to stick with GLM's, it is common for logistic regression to use penalized regression for perfectly separated data. It may also make sense here for Poisson, although there isn't an obvious default for the penalties, nor am I sure about how the standard errors/coverage for coefficients will behave.