通过一些实验计算添加到其他答案。大的标准误差managementD是由小样本量引起的。您得到的标准误差基于近似值，基于对数似然函数近似二次，但事实并非如此。我们可以尝试通过分析获得置信区间，但 Rprofile函数不适用于glm.nb，所以我尝试使用包的解决方法bbmle：

counts <- (c(67, 194, 155, 135, 146, 257, 114, 134, 111, 87, 
               62, 67, 85, 89, 63, 86, 97, 44, 0, 0, 0, 0, 0, 0))
management <- rep(LETTERS[1:4], each = 6)

mydf <- data.frame(counts, management)

model.bbmle <- mle2(counts ~ dnbinom(mu=exp(logmu), size=exp(logtheta)),
                    method= "BFGS", parameters=list(logmu ~ 0 +  
                       management),
                    data=mydf, start=list(logmu=0, 
                      logtheta=2.42),
                    control=list(trace=1)  )

summary(model.bbmle)
Maximum likelihood estimation

Call:
mle2(minuslogl = counts ~ dnbinom(mu = exp(logmu), size = exp(logtheta)), 
    start = list(logmu = 0, logtheta = 2.42), method = "BFGS", 
    data = mydf, parameters = list(logmu ~ 0 + management), control = list(trace = 1))

Coefficients:
                   Estimate Std. Error z value     Pr(z)    
logmu.managementA   5.06898    0.12585 40.2783 < 2.2e-16 ***
logmu.managementB   4.56265    0.12856 35.4897 < 2.2e-16 ***
logmu.managementC   4.34811    0.13017 33.4038 < 2.2e-16 ***
logmu.managementD -11.55880  132.09516 -0.0875    0.9303    
logtheta            2.42213    0.36435  6.6478 2.976e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: 175.9181 

拟合与您的 相当glm.nb，最大对数似然相等，（注意更改的参数化），标准误差较低，但仍然很大！

现在，我们可以尝试 profiling，但是效果不是很好，所以我只给出代码：

prof.4 <- bbmle::profile(model.bbmle, which=4, maxsteps=1000, 
                 alpha=0.005, trace=TRUE)  

 confint(prof.4)

    2.5 %    97.5 % 
       NA -89.53398 
 Warning messages:
 1: In .local(object, parm, level, ...) :
   non-monotonic spline fit to profile (logmu.managementD):  reverting from spline to linear approximation
 2: In regularize.values(x, y, ties, missing(ties), 
      na.rm = na.rm) :
  collapsing to unique 'x' values

返回的间隔（记住对数刻度！）没有意义，我们很快就会明白为什么。我将沿D轴显示对数似然函数的一部分（这与分析不同，因为其他参数保持固定）。这是一些我不完全理解的丑陋代码（由bbmle使用 S4 对象系统引起）：

B <- coef(model.bbmle)

minuslogl_0 <- slot(model.bbmle, "minuslogl")

minuslogl <- function(B) do.call("minuslogl_0", namedrop(as.list(B)))

但是现在我们可以沿着D轴绘制负对似然函数的一部分，其中其他参数保持在 maxlik 估计值：

轴上x是D参数与其 maxlik 值的偏差。可以看到不能设置下限（或者，在原始尺度上，0 是下限），但是可以设置一个尖锐的上限，它将小于标准误差计算所指示的值。使用的代码是

delta <- 10
plot( Vectorize( function(x) minuslogl(B + c(0, 0, 0, x, 0)) ),  
     from=-delta, to=delta, ylab="minusloglik", 
     main="Section of negative loglikelihood function",
     col="red")

我对@dariober 的第一个回答有些异议。

1. 添加1是软糖。

2. 没有实质性的理由不相信样本中记录的零。

3. 最重要的是，模型拟合是合理的，唯一奇怪的是在一种情况下置信区间相当宽。有一些稳健性，因为泊松和负二项式拟合在拟合值上基本相同。（事实上​​，对于这个结构，所有合理的模型，以及一些不太合理的模型，基本上返回组均值作为拟合值。唯一的区别是推断的小字，如果你对此感到不安，你真的需要一个更大的数据集！简单说……）

一张图说明了一切：

为了完整起见，这是我使用的 Stata 代码。自然，在任何现代统计环境中，计算都很简单。

clear 
mat counts = (67,194,155,135,146,257,114,134,111,87,62,67,85,89,63,86,97,44,0,0,0,0,0,0)
set obs 24
gen counts = counts[1, _n]
egen management = seq(), block(6)
label define management 1 A 2 B 3 C 4 D
label val management management

glm counts i.management , family(poisson)
predict poisson
glm counts i.management , family(nbinomial)
predict nbinomial 

* uncomment next if you need to install 
* ssc install stripplot 

gen management1 = management - 0.1
gen management2 = management - 0.2 

stripplot counts , over(management)  vertical stack height(0.3) legend(on order(1 "data" 2 "Poisson fit" 3 "Negative binomial fit")) yla(, ang(h)) ///
addplot(scatter poisson management2, ms(D) || scatter nbinomial management1, ms(T))

编辑对于注入贝叶斯风味的临时方法比仅将所有计数加 1稍微少一些，我使用了 IJ Good 建议的准贝叶斯平滑（对于独立帐户，请参阅本文；本文中的错字修复，pp. 494-495），然后通过泊松 GLM（使用稳健的（三明治-Huber-Eicker-White）标准误差）推动这些调整后的计数。P 值更有意义，同时预测的平均值与任何其他拟合没有太大区别。将有其他可以说是更好的方法来做到这一点。

-------------------------------------------------------------
  management |       mean     Poisson   nbinomial  qs_Poisson
-------------+-----------------------------------------------
           A |      159.0       159.0       159.0       157.7
           B |       95.8        95.8        95.8        95.6
           C |       77.3        77.3        77.3        77.4
           D |        0.0         0.0         0.0         1.5
-------------------------------------------------------------

按照 OP 在评论中的要求，我将在 R 中举一个例子，按照@GordonSmyth 的建议，应用似然比检验 (LRT) 来测试管理组之间的差异。我不确定我是否做对了，所以请检查一下 - 归功于 Gordon，错误是我的。

使用 LRT，我们检查嵌套模型之间的显着差异。要将其应用于这种情况，我们需要将因子扩展management为矩阵（我想glm无论如何都会在内部这样做）。然后我们可以依次去掉每个因素，看看更简单的模型是否与完整的不同：

library(MASS)

counts <- c(67, 194, 155, 135, 146, 257, 114, 134, 111, 87, 
            62, 67, 85, 89, 63, 86, 97, 44, 0, 0, 0, 0, 0, 0)

management <- rep(LETTERS[1:4], each = 6)

design <- model.matrix(~ management)
design
   (Intercept) managementB managementC managementD
1            1           0           0           0
2            1           0           0           0
3            1           0           0           0
4            1           0           0           0
5            1           0           0           0
6            1           0           0           0
7            1           1           0           0
8            1           1           0           0
9            1           1           0           0
10           1           1           0           0
11           1           1           0           0
12           1           1           0           0
13           1           0           1           0
14           1           0           1           0
15           1           0           1           0
16           1           0           1           0
17           1           0           1           0
18           1           0           1           0
19           1           0           0           1
20           1           0           0           1
21           1           0           0           1
22           1           0           0           1
23           1           0           0           1
24           1           0           0           1

拟合完整模型，我们告诉glm.nb忽略截距，因为这已经编码在设计矩阵中。您可能需要使用以下方法检查是否相同glm.nb(counts ~ management)：

fit_full <- glm.nb(counts ~ 0 + design)

现在我们放弃 B 组，拟合简化模型并与完整模型进行比较。这应该等同于评估差异组 A 和组 B 的显着性。我们得到的 p 值为 ~0.01：

design_red <- design[, - which(colnames(design) == 'managementB')]
fit_red <-  glm.nb(counts ~ 0 + design_red)
anova(fit_full, fit_red)

Likelihood ratio tests of Negative Binomial Models

Response: counts
           Model  theta Resid. df    2 x log-lik.   Test    df LR stat. Pr(Chi)
1 0 + design_red  7.668        21          -182.4                              
2     0 + design 11.274        20          -175.9 1 vs 2     1    6.531  0.0106

我们可以对 D 组做同样的事情：

design_red <- design[, - which(colnames(design) == 'managementD')]
fit_red <-  glm.nb(counts ~ 0 + design_red)
anova(fit_full, fit_red)
Likelihood ratio tests of Negative Binomial Models

Response: counts
           Model     theta Resid. df    2 x log-lik.   Test    df LR stat. Pr(Chi)
1 0 + design_red 829670.23        21         -1640.8                              
2     0 + design     11.27        20          -175.9 1 vs 2     1     1465       0

不出所料，A 和 D 之间差异的 p 值接近 0。请注意，简化模型的 theta 参数很大并且会glm.nb发出警告。我不确定如何解释这些，但我想这并不奇怪，因为截距包括带有一串零的大值。

例如，为了测试 B 和 CI 之间的差异，将重新编码整个矩阵以使用 B 而不是 A 作为截距并按上述方式进行 - 我认为有更好的方法。

希望这会有所帮助，我做对了。但是，我仍然认为我添加伪计数的其他解决方案值得考虑。

有关为什么会发生这种情况的解释，请参阅GLM 以获取一个类别中全为零的计数数据。

我知道我可以做一个贝叶斯方法，但在我开始之前，我想知道是否有一种解决方法可以让我使用常客来做这个分析

我认为您可以将所有观察结果加 1，然后得到合理的估计：

model <- glm.nb(counts+1 ~ management)
> summary(model)

Call:
glm.nb(formula = counts + 1 ~ management, init.theta = 
       11.87310622, link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4588  -0.3605   0.0000   0.4648   1.7339  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   5.0752     0.1228  41.330  < 2e-16 ***
managementB  -0.5022     0.1756  -2.860  0.00424 ** 
managementC  -0.7142     0.1768  -4.041 5.33e-05 ***
managementD  -5.0752     0.4425 -11.470  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

    Null deviance: 310.516  on 23  degrees of freedom
Residual deviance:  18.652  on 20  degrees of freedom
AIC: 198.44

Number of Fisher Scoring iterations: 1

              Theta:  11.87 
          Std. Err.:  4.23 

 2 x log-likelihood:  -188.443 

基本原理是零并不是真正可能的结果，它们只是抽样的结果。通过添加 1，您可以重置为更现实的下限。当然，如果加 1 不会使数据集倾斜太多，这是可以接受的，这似乎不是这里的情况。我认为逻辑与贝叶斯方法并没有太大的不同。

我添加了模拟结果来说明将原始计数加 1 的效果。顺便说一句，我认为这也说明了偏差和方差之间的权衡。

我将 A、B 和 C 组的计数保持不变（我认为这在这里不重要）。对于组 DI，模拟来自具有不同平均值和大小的负二项分布的计数。然后我拟合glm.nb模型加和不加 1。为简单起见，我将 D 组设为截距。

这是 D 组估计值的摘要图，请注意，为了便于可视化，我将低于 -10 的估计值重置为 -10：

水平线是真正的手段。每个盒子是 500 个模拟。

• 当真实平均值非常低（大约 < 1）时，校正计数的估计值明显高估了真实平均值（它们是有偏差的），但它们在重复之间非常一致（低方差）。

• 相反，在均值较低的情况下，原始计数的估计值在重复之间非常不稳定（高方差）。重复相同的实验后，您可能会得到截然不同的结果。

• 当真实均值大于 5 时，加 1 的效果开始消失。

我认为由分析师决定如何进行，但凭直觉，我会说在这里添加 1“有意义”。我认为问题不在于未校正估计的数值稳定性，而在于它们是否有意义。

代码：

library(data.table)
library(ggplot2)

counts <- c(67, 194, 155, 135, 146, 257, 114, 134, 111, 87, 62, 
       67, 85, 89, 63, 86, 97, 44)
management <- rep(LETTERS[1:4], c(6, 6, 6, 6))
management <- relevel(as.factor(management), ref= 'D')

seed <- 1
dat <- list()
for(i in 1:500) {
    for(mu in c(0.1, 1, 5, 10)) {
        for(size in c(0.1, 1, 10)) {
            set.seed(seed)
            d <- data.table(
                seed= seed,
                mu= mu,
                size= size,
                counts= c(counts, rnbinom(n= sum(
                 management == 'D'), mu= mu, size= size)),
                management= management
            )
            dat[[length(dat) + 1]] <- d
            seed <- seed + 1
        }
    }
}
dat <- rbindlist(dat)

estimates <- dat[, list(raw= coef(glm.nb(counts ~   
    management))[1], corrected= coef(glm.nb(counts + 1 ~ 
    management))[1]), by= list(seed, mu, size)]
estimates <- melt(estimates, id.vars= c('seed', 'mu', 'size'), variable.name= 'method', value.name= 'estimate')

estimates[, label := sprintf('True mean count= %s', mu)]
estimates[, label := factor(label, sprintf('True mean count= %s', unique(sort(mu))))]
gg <- ggplot(data= estimates, aes(x= as.factor(size), y= 
      ifelse(estimate < -10, -10, estimate), colour= method)) +
        geom_hline(data= unique(estimates[,list(mu, label)]), 
      aes(yintercept= log(mu)), colour= 'grey30', 
      linetype= 'dashed') +
        geom_boxplot() +
        xlab('Size (dispersion parameter)') +
        ylab('Estimate (capped to -10)') +
        theme_light() +
        theme(strip.text= element_text(colour= 'black', 
        size= 12)) +
        facet_wrap(~ label, scales= 'free_y')

编辑：我对此解决方案持怀疑态度，对此我感到有些惊讶，因为在计数数据中添加 1 并非闻所未闻。

例如，来自测序技术的基因表达数据就是计数​​数据。一种非常受人尊敬的差异表达分析方法是 limma-voom（Gordon Smyth 是合著者，如果我引用错误，我深表歉意）。从论文：

计数从零偏移 0.5 以避免取零的对数，并减少低表达基因的 log-cpm 可变性。

这正是我在这里提出的。

基因表达永远不会真正为零，通常情况下，计数 0 和计数 1 之间的差异在生物学上是无关紧要的。OP的情况是否相同？我不知道，但我怀疑它是。D 组有 0 条斑马，还是只有 1 条或 2 条在 6 次观察中未被检测到，这真的很重要吗？如果答案是否定的，那么加 1 比产生 -Inf 系数更明智。