Does glm.nb (in R) give inflated p-values when tested under the null? I'm trying to do a simple test of glm.nb in R.
I simulate outcomes from a negative binomial using rnegbin.  I have a 0/1 case/control variable, and I'd like to test it for significance, but I'm finding that there are too many significant p values when I test null data:
pvalues = c()
for (blah in 1:20000) {
  outcomes = rnegbin(20,mu=30,theta=5)
  casecontrol =  c(rep(0,10),rep(1,10))
  model = glm.nb(outcomes~casecontrol,maxit=1000)
  pvalues = c(pvalues,summary(model)$coefficients[2,4])
}
hist(pvalues,40)


What am I doing wrong here?  I have negative binomial data and I'm testing it with negative binomial regression -- maybe there's something fundamental I'm not understanding?  Would be hugely appreciative if someone could point me in the right direction.
 A: The usual p-values based on maximum likelihood methods typically involve t statistics constructed by dividing the estimated parameters by the (Wald) standard errors, and then comparing the results to the standard normal distribution.  There are two sources of error here: The first is that, like the usual regression model, there is variability in the standard error that makes the t distribution more appropriate than the z distribution; and second error is that the distribution of the estimate is only approximately normal. Both of these problems lessen with larger sample sizes.  Here is some modified code that addresses both the t vs z and the sample sizes issues.
pvalues = c()
tvalues = c()
ndiv2 = 10
for (blah in 1:20000) {
outcomes = rnegbin(2*ndiv2,mu=30,theta=5)
casecontrol =  c(rep(0,ndiv2),rep(1,ndiv2))
model = glm.nb(outcomes~casecontrol,maxit=1000)
pvalues = c(pvalues,summary(model)$coefficients[2,4])
tvalues = c(tvalues,summary(model)$coefficients[2,3])
}
pvalues1 = 2*(1 - pt(abs(tvalues), 2*(ndiv2-1)))
hist(pvalues,40)
hist(pvalues1, 40)
mean(pvalues <=0.05)  # should be close to .05
mean(pvalues  <= 0.01)  # should be close to .01
mean(pvalues  <=0.005)  # should be close to .005
mean(pvalues <= 0.001)  # should be close to .001
mean(pvalues1 <=0.05)  # should be close to .05
mean(pvalues1  <= 0.01)  # should be close to .01
mean(pvalues1  <=0.005)  # should be close to .005
mean(pvalues1 <= 0.001)  # should be close to .001

Even with your small sample size of 20 (and to be clear, that is the issue with your results), the t-based results look a lot better.  And if you increase the sample size from 20 to 200 (by changing ndiv2 to 100), the results look even better.
