glmmTMB拟合Poisson混合模型遇NaN与收敛问题求助

我用R语言的glmmTMB包拟合Poisson分布的广义线性混合模型，代码如下：

mod_po <- glmmTMB(richness ~ x_dif + x_dif2 + (1|papers/watersheds), data = df_glm, family = poisson)

其中x_dif是河流尺寸的线性变量，x_dif2是它的二次项，已经对变量做了中心化处理，但二者相关度仍高达0.85。由于研究假设要求，不能移除二次项。

拟合模型时出现警告：

In sqrt(diag(vcovs)) : NaNs produced

模型摘要里x_dif2的标准误等为NaN；计算R²时收到警告：

Model convergence problem; non-positive-definite Hessian matrix

模型摘要信息 #

> summary(mod_po)
 Family: poisson  ( log )
Formula:          
riqueza_larguras_unl ~ x_dif + x_dif2 + (1 | nomes_artigos_largura_unl/bacias)
Data: df_glm
     AIC      BIC   logLik deviance df.resid 
  2927.8   2946.6  -1458.9   2917.8      314 
Random effects:
Conditional model:
 Groups                           Name        Variance Std.Dev.
 bacias:nomes_artigos_largura_unl (Intercept) 0.09382  0.3063  
 nomes_artigos_largura_unl        (Intercept) 0.09382  0.3063  
Number of obs: 319, groups:  
bacias:nomes_artigos_largura_unl, 24; nomes_artigos_largura_unl, 24
Conditional model:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.087e+00  8.982e-02   34.37  < 2e-16 ***
x_dif       -7.929e-03  1.188e-03   -6.67 2.53e-11 ***
x_dif2       2.624e-05        NaN     NaN      NaN    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In sqrt(diag(vcovs)) : NaNs produced

诊断信息 #

predictors with unusually large or small standard deviations (|log10(sd)|>2120.57):
  x_dif2 
2120.565 
Predictor variables with very narrow or wide ranges generally
give rise to parameters with very large or small magnitudes,
which can sometimes exacerbate numerical instability, and may
also be appear (incorrectly) to be indicating a poorly defined
optimum (i.e., a non-positive definite Hessian
Unusually large Z-statistics (|x|>5):
(Intercept)       x_dif 
  34.368298   -6.671379 
Large Z-statistics (estimate/std err) suggest a *possible*
failure of the Wald approximation - often also associated with
parameters that are at or near the edge of their range (e.g.
random-effects standard deviations approaching 0).
(Alternately, they may simply represent very well-estimated
parameters; intercepts of non-centered models may fall in this
category.) While the Wald p-values and standard errors listed
in summary() may be unreliable, profile confidence intervals
(see ?confint.glmmTMB) and likelihood ratio test p-values
derived by comparing models (e.g. ?drop1) are probably still
OK.  (Note that the LRT is conservative when the null value is
on the boundary, e.g. a variance or zero-inflation value of 0
(Self and Liang 1987; Stram and Lee 1994; Goldman and Whelan
2000); in simple cases these p-values are approximately twice
as large as they should be.)

计算R²时的警告信息 #

1: In r.squaredGLMM.glmmTMB(mod_po) :
  the effects of zero-inflation and dispersion model are ignored
2: the null model is correct only if all variables used by the original model remain unchanged. 
3: In fitTMB(TMBStruc) :
Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting').

请问该问题是否由变量相关性导致？如何在保留二次项的前提下解决？

这个问题的核心原因是线性项与二次项的高度共线性，再加上x_dif2的数值尺度异常大（诊断信息显示其log10标准差远超阈值），双重因素引发了模型的数值不稳定，最终导致Hessian矩阵非正定、标准误出现NaN。以下是保留二次项的可行解决方法：

1. 标准化自变量（替代仅中心化） #

中心化只能消除截距的解释偏差，但无法缓解共线性。标准化（将变量转换为均值为0、标准差为1的z-score）可以大幅缩小二次项的数值尺度，同时降低线性项与二次项的相关程度。操作代码：

# 标准化线性项
df_glm$x_dif_z <- scale(df_glm$x_dif)
# 基于标准化后的变量计算二次项，避免手动计算引入的尺度误差
df_glm$x_dif2_z <- df_glm$x_dif_z^2

# 重新拟合模型
mod_po_z <- glmmTMB(richness ~ x_dif_z + x_dif2_z + (1|papers/watersheds), 
                    data = df_glm, family = poisson)

2. 使用正交多项式替代原始二次项 #

正交多项式可以让线性项和二次项完全不相关，从根源上解决共线性问题。R中用poly()函数生成正交多项式即可：

# 生成2阶正交多项式（raw=FALSE为正交模式，默认开启）
poly_terms <- poly(df_glm$x_dif, degree = 2, raw = FALSE)
df_glm$poly1 <- poly_terms[,1]  # 对应线性效应
df_glm$poly2 <- poly_terms[,2]  # 对应二次效应

# 拟合模型
mod_po_poly <- glmmTMB(richness ~ poly1 + poly2 + (1|papers/watersheds), 
                       data = df_glm, family = poisson)

注意：正交多项式的系数解释和原始二次项不同，但可以通过转换还原回原始尺度的效应，完全不影响“是否存在二次效应”的假设检验。

3. 调整模型优化参数 #

如果尺度调整后仍有收敛问题，可以尝试更换优化器或增加迭代次数：

# 使用BFGS优化器并增加迭代次数
mod_po_optim <- glmmTMB(richness ~ x_dif_z + x_dif2_z + (1|papers/watersheds), 
                        data = df_glm, family = poisson,
                        control = glmmTMBControl(optimizer = optim, 
                                                 optArgs = list(method = "BFGS", maxit = 1000)))

# 或者调整nlminb优化器的收敛阈值
mod_po_tight <- glmmTMB(richness ~ x_dif_z + x_dif2_z + (1|papers/watersheds),
                        data = df_glm, family = poisson,
                        control = glmmTMBControl(optArgs = list(eval.max = 1000, iter.max = 1000)))

4. 改用稳健的统计推断方法 #

如果模型仍存在标准误NaN的情况，不要依赖模型摘要中的Wald检验（z值、p值），改用以下方法：

• 似然比检验：用drop1()比较包含/不包含二次项的模型，检验二次效应的显著性：

  drop1(mod_po_z, test = "Chisq")
  

• 轮廓置信区间：用confint()生成二次项系数的置信区间，替代标准误：

  confint(mod_po_z, parm = "x_dif2_z")
  

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