Instead of overdispersed (or quasi-)poisson regression you can use the NB1 distribution, which has the same linear variance function as ODP and a full-fledged likelihood function instead of the quasilikelihood of ODP. NB1 is implemented in the gamlss
package as family=NBII
, whereas regular Negative Binomial can be called through family=NBI
. All credit for this part of the answer goes to @Achim Zeileis for helping me with a similar question here: Why is the Quasipoisson in glm not treated as a special case of Negative Binomial? , see his post for more info regarding NB/NB1 (and the confusing naming conventions).
Regarding ANOVA, I have not been able to find a built-in method for gamlss
objects, but it is not hard to write your own implementation of the chi-squared test statistic. An example:
set.seed(123)
data = rNBII(100,mu = 6,sigma=0.5) #generate NB1 data with mean mu=5 and variance (1+sigma)*mu = 9
h0 = gamlss(data~1,family=PO) #null model: poisson
h1 = gamlss(data~1,family=NBII) #alternative model: NB1/ODP
df = h1$df.fit - h0$df.fit
deviance = as.numeric(-2*logLik(h0) + 2*logLik(h1))
p.value = pchisq(deviance,df,lower.tail=F)
> p.value
[1] 0.01429169 #reject the null model at > 95% confidence