如何为R中MASS包的glm.nb函数选择正确的init.theta初始值

1. init.theta的具体含义 #

init.theta是负二项分布离散参数θ的迭代初始值。glm.nb()拟合模型时，需要通过极大似然估计迭代优化得到最终的θ值，这个迭代过程的起点就是init.theta的传入值。
默认情况下，函数会先拟合泊松GLM，再基于残差计算初始θ值，但如果你的数据离散程度和泊松假设偏差过大，默认计算的初始值会离真实最优θ很远，就会出现迭代上限触达的报错，同时输出的θ、标准误、p值都会失真。

2. 合理init.theta的确定方法 #

• 预计算估计法
  先通过泊松模型粗略估算数据的离散程度，再推导初始θ，避免初始值偏离过大：

  # 拟合基准泊松模型
  poi_mod <- glm(measurement ~ year, data = df, family = poisson(log))
  # 计算离散系数：皮尔逊卡方/自由度，离散系数>1说明数据过离散
  disp_ratio <- sum(residuals(poi_mod, type = "pearson")^2) / df.residual(poi_mod)
  # 推导初始θ：θ和离散系数成反比
  init_theta <- 1 / disp_ratio
  # 代入拟合负二项模型
  nb_mod <- glm.nb(measurement ~ year, data = df, init.theta = init_theta, link = log)
  

• 网格试错法
  如果预计算的初始值依然无法收敛，可以在合理范围（通常是0.1~10，过离散越严重初始值可以设得越小）内尝试多个初始值，选择无收敛警告、AIC最低的模型结果即可。手动设置的init.theta = 1是通用的常用初始值，只要拟合后无报错警告，结果就是可信的。
• 调整迭代上限法
  也可以不手动设置初始值，直接提高函数的迭代次数上限，给算法足够的迭代空间找到最优θ：

  nb_mod <- glm.nb(measurement ~ year, data = df, link = log, control = glm.control(maxit = 100))
  

3. 结果合理性校验 #

最终拟合的模型满足以下两个条件，即可确认θ取值合理：

• 模型运行无迭代上限、收敛失败类警告
• 模型残差无明显的异常趋势，离散程度符合预期

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