如何计算负二项GLM所需theta值并解决glm.nb运行报错问题

问题原因 #

从输出可以看到，模型已经完成了初步拟合，估计得到的theta值为2055606，后续报错是summary函数格式化输出时触发的，两个警告的核心原因是theta估计的默认迭代次数不足，且theta估计值过大导致数值计算出现NaN。
theta值越大代表负二项分布越接近泊松分布，你当前的theta值已经超过200万，说明数据的过度离散程度非常弱，接近泊松分布的适用条件。

解决方案 #

• 方案1：直接从模型对象提取theta值
  不需要等待summary函数输出完成，拟合完成后直接从模型对象中提取参数即可，代码如下：

library(MASS)
m1 <- glm.nb(Total ~ Treatment, data = twohour)
# 直接提取theta
theta_value <- m1$theta
print(theta_value)

• 方案2：提高迭代上限解决警告问题
  glm.nb默认的迭代上限为10次，无法满足大theta的估计需求，手动调大迭代次数即可：

m1 <- glm.nb(Total ~ Treatment, data = twohour, control = glm.control(maxit = 100))
summary(m1)

• 方案3：手动指定初始theta值
  从你已有的运行结果中已经得到初始theta估计值2055606，直接传入作为初始值可以大幅降低迭代需求：

m1 <- glm.nb(Total ~ Treatment, data = twohour, init.theta = 2055606, control = glm.control(maxit = 100))
summary(m1)

• 方案4：替代模型选择
  如果验证后发现数据的过度离散程度极低（可通过泊松模型的离散参数验证，离散参数=残差偏差/自由度，若接近1则离散程度弱），可以直接改用准泊松GLM，不需要估计theta参数即可控制过度离散，结果和当前负二项模型几乎一致：

# 准泊松模型
quasi_poi_m <- glm(Total ~ Treatment, data = twohour, family = quasipoisson)
# 查看离散参数
quasi_poi_m$dispersion

内容的提问来源于stack exchange，提问作者BC2504