零膨胀负二项模型拟合出现NaN警告的原因及解决方法

我尝试使用零膨胀负二项模型（zero inflated negative binomial model）拟合投诉数据及相关预测变量，普通负二项模型运行正常，但零膨胀版本的模型摘要中，visits和hours两个预测变量的标准误（SE）与p值出现NaN值，相关警告及模型输出如下：

Warning: NaNs produced
Call:
zeroinfl(formula = complaints ~ visits + gender + revenue + residency + 
    hours, data = comp, dist = "negbin")

Pearson residuals:
     Min       1Q   Median       3Q      Max 
-1.07975 -0.65173 -0.24077  0.08785  3.31293 

Count model coefficients (negbin with log link):
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.7050283  1.4261221   2.598  0.00938 ** 
visits       0.0022846        NaN     NaN      NaN    
genderM      0.5783319  0.3036332   1.905  0.05682 .  
revenue     -0.0174033  0.0043892  -3.965 7.34e-05 ***
residencyY  -0.1785727  0.3117556  -0.573  0.56678    
hours       -0.0029508  0.0007413  -3.980 6.88e-05 ***
Log(theta)   0.5487658        NaN     NaN      NaN    

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 22.08077    4.43783   4.976 6.51e-07 ***
visits       0.01180        NaN     NaN      NaN    
genderM      0.26974    1.55676   0.173  0.86244    
revenue     -0.11795    0.01141 -10.339  < 2e-16 ***
residencyY   4.56503    1.45564   3.136  0.00171 ** 
hours       -0.01931        NaN     NaN      NaN    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 1.7311 
Number of iterations in BFGS optimization: 52 
Log-likelihood: -127.2 on 13 Df 

以下是部分0值样本数据（变量顺序为visits、complaints、residency、gender、revenue、hours）：

201   3374          0         Y      M 255.6324 1657.8924
211   1555          0         Y      M 342.8607 1212.1892
221   2185          0         Y      M 294.5911  932.0161
231    972          0         Y      M 311.6149 1070.7129
241   1566          0         Y      M 260.6672 1707.1469
251   2574          0         Y      M 259.4384 1449.8738
261   2390          0         N      F 219.9528 1759.5416
271   2445          0         N      F 243.4169 1156.7693
281   3077          0         N      F 224.2506 1114.2919
291   2908          0         N      F 257.8778 1065.6133
301   3745          0         N      F 246.1194 2269.0726

问题 #

请问该问题的成因是什么？有哪些解决方法？

补充问题 #

我尝试将visits变量缩放（除以10，因数值在千级），NaN警告消失，但此举是否需要对模型或数据做进一步调整以保证一致性？

问题成因 #

• 数值尺度差异引发优化不稳定：visits属于千级变量，和其他变量尺度差距过大，BFGS优化算法在更新参数时难以平衡不同尺度变量的步长，导致用于计算标准误的海森矩阵出现奇异或不可逆情况，最终生成NaN。
• 准完全分离或多重共线性：虽然未发现完全分离，但可能存在准完全分离（某变量对0/非0分组区分度极高但不绝对），或者visits与hours等变量存在高度共线性，使得模型无法估计可靠的标准误。
• 模型参数冗余：零膨胀模型同时拟合计数和零膨胀两个部分，若visits/hours在两个部分的作用高度相关，可能导致参数估计自由度不足，引发数值问题。

解决方法 #

• 变量标准化/缩放：如你已尝试的，对大尺度变量进行缩放（除以10、标准化为均值0方差1等），平衡变量间尺度，帮助优化算法稳定收敛。
• 检查并处理共线性：计算变量间的方差膨胀因子（VIF），若VIF>10说明存在严重共线性，可移除高度相关变量，或使用Lasso等正则化方法约束参数。
• 简化模型结构：先尝试拟合零膨胀泊松模型，再逐步过渡到负二项；或测试仅在计数部分/零膨胀部分包含visits/hours，观察参数估计是否稳定。
• 更换优化算法：默认BFGS对部分数值问题敏感，可通过zeroinfl的control参数指定其他优化器，例如：

  zeroinfl(formula = complaints ~ visits + gender + revenue + residency + hours, 
           data = comp, dist = "negbin", 
           control = list(method = "Nelder-Mead"))
  

• 处理异常值与分布偏态：检查visits和hours的分布，若存在极端偏态或异常值，可对极端值进行截断或对数转换，减少对优化过程的干扰。

关于变量缩放的后续处理 #

• 无需调整模型结构：变量缩放仅改变参数尺度，不影响模型拟合效果和预测能力。例如visits除以10后，参数估计值变为原来的10倍，解释时只需对应缩放后的变量（每增加10次访问对应原参数的效应）。
• 保持解释一致性：报告结果时需明确说明变量缩放方式，避免解释错误。比如原visits每增加1单位的效应，对应缩放后变量每增加10单位的效应。
• 建议标准化所有连续变量：为统一解释逻辑，可对所有连续预测变量（visits、revenue、hours）进行标准化，此时参数可表示为变量每变化1个标准差对应的效应，便于比较不同变量的影响大小。

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