y <- c(5, 12, 25, 6, 10, 8, 4, 3, 0, 0, 1, 0, 0, 0, 0, 0)
fac <- rep(LETTERS[1:4], each = 4)
library(MASS)
model <- glm.nb(y ~ fac)
summary(model)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    2.4849     0.2749   9.041  < 2e-16 ***
facB          -0.6523     0.4126  -1.581 0.113899
facC          -3.8712     1.0631  -3.641 0.000271 ***
facD         -21.7875  4713.3084  -0.005 0.996312
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(4.5692) family taken to be 1)