glm.nb()大数据集报错：无法找到有效系数，如何设置初始值？

问题背景 #

使用glm包的glm.nb()拟合负二项回归模型：

m_nb= glm.nb(Error_Count ~ TotalWL + Auto_frac + PHONE + JUSTIF_weight + MESSAGE_OTHER_count + Hour +
I(Auto_frac^2)+I(TotalWL^2), data = df)

10000条数据集可正常运行，但60000条大数据集时报错：

Error: no valid set of coefficients has been found: please supply starting values

手动指定初始值后仍报错：

Error: cannot find valid starting values: please specify some

已知fixest包的fenebin()可成功拟合该模型，但因需要predict()输出的标准误，必须使用glm包。

解决方案：从fixest模型提取初始值 #

既然fenebin()能收敛，直接用它的拟合结果作为glm.nb()的初始值，步骤如下：

1. 用fixest拟合模型并提取参数

# 加载fixest包
library(fixest)
# 拟合负二项模型
m_fix <- fenebin(Error_Count ~ TotalWL + Auto_frac + PHONE + JUSTIF_weight + MESSAGE_OTHER_count + Hour + I(Auto_frac^2) + I(TotalWL^2), data = df)
# 提取回归系数（顺序与公式变量严格对应）
coef_vals <- coef(m_fix)
# 提取负二项模型的分散参数theta
theta_val <- m_fix$theta
# 组合成glm.nb需要的初始值向量：先回归系数，最后是theta
start_vals <- c(coef_vals, theta_val)

2. 用提取的初始值运行glm.nb

library(MASS)
# 拟合模型，传入初始值，同时可调整迭代参数提升收敛概率
m_nb <- glm.nb(Error_Count ~ TotalWL + Auto_frac + PHONE + JUSTIF_weight + MESSAGE_OTHER_count + Hour + I(Auto_frac^2) + I(TotalWL^2), 
               data = df, 
               start = start_vals,
               control = glm.control(maxit = 1000, epsilon = 1e-6))

额外优化建议 #

• 中心化变量减少共线性：原模型包含变量及其平方项，易出现多重共线性，可先对变量中心化后再构造平方项，降低共线性帮助收敛：

# 中心化变量
df$TotalWL_c <- df$TotalWL - mean(df$TotalWL, na.rm = TRUE)
df$Auto_frac_c <- df$Auto_frac - mean(df$Auto_frac, na.rm = TRUE)
# 用中心化后的变量重新拟合模型
m_nb <- glm.nb(Error_Count ~ TotalWL_c + Auto_frac_c + PHONE + JUSTIF_weight + MESSAGE_OTHER_count + Hour + I(Auto_frac_c^2) + I(TotalWL_c^2), 
               data = df, 
               start = start_vals,
               control = glm.control(maxit = 1000, epsilon = 1e-6))

• 检查极端值：大数据集中的极端值可能干扰收敛，可先对连续变量（如TotalWL、JUSTIF_weight）做极值截断处理，比如截去上下1%的极端值。

内容的提问来源于stack exchange，提问作者georgia-max