【问题标题】：Scaling predictors in lme4 glmer doesn't resolve eigenvalue warnings; neither does alternative optimizationlme4 glmer 中的缩放预测变量不能解决特征值警告；替代优化也没有
【发布时间】：2018-12-22 04:12:57
【问题描述】：

我正在使用 R 中的lme4 的glmer 函数分析数据（包括在下面）。我正在构建的模型包含一个泊松分布响应变量 (obs)、一个随机因子 (area)、一个连续偏移 (duration)、五个连续固定效应 (can_perc、can_n、 time、temp、cloud_cover）和一个二项式固定效应因子（burnt）。在拟合模型之前，我检查了共线性并删除了所有共线性变量。

初始模型为：

q1 = glmer(obs ~ can_perc + can_n  + time * temp + 
           cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), 
           data=dat, family=poisson, na.action = na.fail)

（注意：我需要将na.action 指定为“na.fail”，因为我想稍后将dredge() 模型指定为这是必需的。）

运行模型会给出以下警告：

“Hessian 在数值上是奇异的：参数不是唯一确定的”

在模型的类似变体中，我也收到了警告：

"在 checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : 模型几乎无法识别：大特征值比 - 重新调整变量？”

根据我对https://rdrr.io/cran/lme4/man/troubleshooting.html 和其他地方的建议的有限理解，这两个警告都反映了一个类似的问题，即 Hessian（逆曲率矩阵）具有较大的特征值，表明（在数值公差范围内）表面完全在某个方向平坦。根据警告和链接中的建议，我使用scale() 重新调整了所有连续预测变量。我还缩放了偏移变量（我尝试了缩放和不缩放这个变量）。具有比例预测变量的模型在这里：

q2 = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + 
           s.cloud_cover + factor(burnt) + (1|area) +
           offset(dat$s.duration), 
           data=dat, family=poisson, na.action = na.fail)

但是我还没有逃脱特征值！缩放模型给出了两个警告：

“无法评估缩放渐变”
“模型收敛失败：具有1个负特征值的退化Hessian”

我在网上搜索了很多，除了检查模型没有被错误指定之外，找不到其他案例/解决方案来处理一旦预测变量被缩放后如何处理特征值问题。

尝试解决警告/改进优化：

基于这些页面/文档： https://cran.r-project.org/web/packages/lme4/lme4.pdf

https://stats.stackexchange.com/questions/164457/r-glmer-warnings-model-fails-to-converge-model-is-nearly-unidentifiable

https://rdrr.io/cran/lme4/man/isSingular.html

https://stats.stackexchange.com/questions/242109/model-failed-to-converge-warning-in-lmer

和其他人，

我有：

检查模型规格和数据是否有错误（我看不到 - 我错过了什么吗？）
用 is_singular(x, tol = 1e-05) 检查奇异性（不知何故，这个函数调用从 isSingular() 演变为当前形式？）：所有模型都给出 FALSE。
用converge_ok(q2, tolerance = 0.001)检查收敛度量：所有模型都给出FALSE，除非我大幅增加容差；但是它们的收敛度量确实存在很大差异。
尝试了不同的优化器/模型估计方法如下：
- a) glmerControl(optimizer = "bobyqa") and glmerControl(optimizer ="Nelder_Mead")
- b) glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb'))
- c) bobyqa、Nelder_Mead、optimx.nlminb、optimx.L-BFGS-B、nloptwrap.NLOPT_LN_NELDERMEAD、nloptwrap.NLOPT_LN_BOBYQA 和 nmkbw 优化器，使用 optimx 包中的 all_fit() 函数。

代码如下：

# singularity and convergence for first two models:
is_singular(s1, tol = 1e-05) # FALSE (a good thing?)
converge_ok(s1, tol = 1e-05) # FALSE (a bad thing?) 0.0259109730912352

is_singular(s2, tol = 1e-05) # FALSE (a good thing?)
converge_ok(s2, tol = 1e-05) # FALSE (a bad thing?) 0.0023434329028163
# I looked at singularity and converge measures for the others below, but omitted for brevity.

# Alternate optimisations for q1:
q1.bobyqa = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
# Warning 1: unable to evaluate scaled gradient
# Warning 2: Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

q1.neldermead = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) +  (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail,  glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
# Warning: unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined

q1.nlminb = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
# Warning: Parameters or bounds appear to have different scalings. This can cause poor performance in optimization. 
# It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.convergence code 9999 from optimxError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,  :   (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

all_fit(q1)

# Alternate optimisations for q2:
q2.bobyqa = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) +  (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
# Warning 1: unable to evaluate scaled gradient
# Warning 2: Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

q2.neldermead = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) + (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
# Warning: unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined

q2.nlminb = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) +  (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, control = glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
# Warning: Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?

all_fit(q2)

上述代码的输出，对于未缩放的模型 (q1)：

is_singular(s1, tol = 1e-05) # FALSE (a good thing?)
[1] FALSE
converge_ok(s1, tol = 1e-05) # FALSE (a bad thing?) 0.0259109730912352
0.0259109730912352 
             FALSE 
is_singular(s2, tol = 1e-05) # FALSE (a good thing?)
[1] FALSE
alternate optimisations for original model:
q1.bobyqa = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues

alternate optimisations for original model:
q1.bobyqa = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
    unable to evaluate scaled gradientModel failed to converge: degenerate      Hessian with 1 negative eigenvalues
    q1.neldermead = glmer(obs ~ can_perc + can_n  + time * temp + cloud_cover + factor(burnt) + (1|area) + offset(dat$duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
    unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined

all_fit(q1)
bobyqa. : unable to evaluate scaled gradientModel failed to converge:     degenerate  Hessian with 1 negative eigenvalues[OK]
Nelder_Mead. : unable to evaluate scaled gradient Hessian is numerically singular: parameters are not uniquely determined[OK]
optimx.nlminb : Parameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA,     NEWUOA.convergence code 9999 from optimxParameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA,     NEWUOA.convergence code 9999 from optimx[ERROR]
optimx.L-BFGS-B : Parameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.convergence code 9999 from optimxParameters or bounds appear to have different scalings.
This can cause poor performance in optimization. 
It is important for derivative free methods like BOBYQA, UOBYQA,     NEWUOA.convergence code 9999 from optimx[ERROR]
nloptwrap.NLOPT_LN_NELDERMEAD : [ERROR]
nloptwrap.NLOPT_LN_BOBYQA : [ERROR]
nmkbw. : [ERROR]

$`bobyqa.`
    Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
     Family: poisson  ( log )
    Formula: obs ~ can_perc + can_n + time * temp + cloud_cover + factor(burnt) +  (1 | area) + offset(dat$duration)
       Data: dat
      AIC       BIC    logLik  deviance  df.resid 
     311.0473  330.3356 -146.5237  293.0473        54 
    Random effects:
     Groups Name        Std.Dev.
     area   (Intercept) 1.992   
    Number of obs: 63, groups:  area, 8
    Fixed Effects:
(Intercept)  can_perc  can_n   time     temp  
-67.4998    -1.3180    0.0239    4.8025    1.7793  
cloud_cover  factor(burnt)unburnt             time:temp  
             -0.3813               18.5676               -0.1748  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$Nelder_Mead.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
     Family: poisson  ( log )
    Formula: obs ~ can_perc + can_n + time * temp + cloud_cover + factor(burnt) +      (1 | area) + offset(dat$duration)
       Data: dat
      AIC       BIC    logLik  deviance  df.resid 
     311.0473  330.3356 -146.5237  293.0473        54 
    Random effects:
     Groups Name        Std.Dev.
     area   (Intercept) 1.992   
    Number of obs: 63, groups:  area, 8
    Fixed Effects:
             (Intercept)              
can_perc      can_n       time      temp  
-67.48057    -1.31791    0.02389    4.80463    1.78012  
cloud_cover  factor(burnt)unburnt    time:temp  
-0.38118    18.52637      -0.17483  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$optimx.nlminb
<std::runtime_error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat,    compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate>

$`optimx.L-BFGS-B`
<std::runtime_error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate>

$nloptwrap.NLOPT_LN_NELDERMEAD
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nloptwrap.NLOPT_LN_BOBYQA
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nmkbw.
<std::runtime_error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate>

上述代码的输出，对于缩放模型 (q2)：

alternate optimisations for q2:
q2.bobyqa = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) + (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5)))
Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?
q2.neldermead = glmer(obs ~ s.can_perc + s.can_n  + s.time * s.temp + s.cloud_cover + factor(burnt) + (1|area) + offset(dat$s.duration), data=dat, family=poisson, na.action = na.fail, glmerControl(optimizer ="Nelder_Mead", optCtrl = list(maxfun = 2e5)))
unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues

all_fit(q2)
bobyqa. : Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?[OK]
Nelder_Mead. : unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues[OK]
optimx.nlminb : Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?[OK]
optimx.L-BFGS-B : unable to evaluate scaled gradientModel failed to converge: degenerate  Hessian with 1 negative eigenvalues[OK]
nloptwrap.NLOPT_LN_NELDERMEAD : [ERROR]
nloptwrap.NLOPT_LN_BOBYQA : [ERROR]
nmkbw. : [ERROR]
$`bobyqa.`
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8412  336.1294 -149.4206  298.8412        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.523   
Number of obs: 63, groups:  area, 8
Fixed Effects:
(Intercept)    s.can_perc    s.can_n    s.time    s.temp  
-18.19816    -0.22152    0.45839    0.05239    -0.24983  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19691              17.92390              -0.13948  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 

$Nelder_Mead.
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8408  336.1290 -149.4204  298.8408        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.524   
Number of obs: 63, groups:  area, 8
Fixed Effects:
         (Intercept)            s.can_perc               s.can_n                s.time                s.temp  
           -19.29632              -0.22153               0.45840               0.05241              -0.24990  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19692              19.02091              -0.13949  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$optimx.nlminb
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8412  336.1294 -149.4206  298.8412        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.523   
Number of obs: 63, groups:  area, 8
Fixed Effects:
         (Intercept)            s.can_perc               s.can_n                s.time                s.temp  
           -18.23626              -0.22152               0.45839               0.05239              -0.24983  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19691              17.96199              -0.13948  
convergence code 0; 1 optimizer warnings; 0 lme4 warnings 

$`optimx.L-BFGS-B`
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: n_shreiberi ~ s.can_perc + s.can_n + s.time * s.temp + s.cloud_cover +  
    factor(burnt) + (1 | area) + offset(dat$s.duration)
   Data: dat
      AIC       BIC    logLik  deviance  df.resid 
 316.8412  336.1294 -149.4206  298.8412        54 
Random effects:
 Groups Name        Std.Dev.
 area   (Intercept) 2.524   
Number of obs: 63, groups:  area, 8
Fixed Effects:
         (Intercept)            s.can_perc               s.can_n                s.time                s.temp  
           -18.23581              -0.22155               0.45841               0.05242              -0.24997  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19694              17.96246              -0.13943  
convergence code 0; 2 optimizer warnings; 0 lme4 warnings 

$nloptwrap.NLOPT_LN_NELDERMEAD
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nloptwrap.NLOPT_LN_BOBYQA
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

$nmkbw.
<simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite>

数据：

数据集可在此链接获得： https://www.dropbox.com/s/ud50uatztjq4bh9/20181217%20Surveys%20simplified%20data%20for%20stackX.xlsx?dl=0

结论与要求：

在我看来，这些替代的优化方法都没有成功；事实上，其中一些似乎引发了其他警告/错误，这将使我陷入另一个兔子洞。

谁能告诉我如何在拟合这些模型方面取得进展？我的目的不是让这些成为最终模型，而是挖掘它们，然后从不同的替代子集模型中选择最佳/顶级模型。

【问题讨论】：

标签： r optimization lme4

【解决方案1】：

tl;dr 这看起来像是完全分离的情况；在你的“烧伤”状态下，你根本没有积极的结果。您不必一定担心这一点 - AIC 比较应该仍然相当可靠 - 但您可能想在继续之前了解正在发生的事情。这个问题（和补救措施）在GLMM FAQ 的相关部分进行了讨论（CrossValidated 上有各种相关的问题/答案）。

我怎么知道？以下是系数：

  (Intercept)       s.can_perc               s.can_n                s.time                s.temp  
   -19.29632          -0.22153               0.45840               0.05241              -0.24990  
       s.cloud_cover  factor(burnt)unburnt         s.time:s.temp  
            -0.19692              19.02091              -0.13949

任何时候（二项式或泊松）GLM 中的系数（绝对值）大于 8-10，您都必须担心（除非您正在查看以非常大的方式测量的数值协变量的系数单位，例如，如果您以千兆吨为单位查看后院的碳含量）。这意味着预测变量中一个单位的变化会导致对数赔率（例如）10 个单位的变化（对于二项式/logit-link 模型），例如从 0.006 (plogis(-5)) 到 0.994 (plogis(5)) 的概率。在您的情况下，截距为 -19.29，因此在所有预测变量的值为零时在燃烧条件下，您得到的概率为 4.2e-9。另一个巨大的系数是unburnt (19.02)，因此在未燃烧（未燃烧？）条件下所有预测变量的值为零时，您会得到 plogis(-19.29+19.02) = 0.43。

【讨论】：

非常感谢 Ben Bolker！