lme4 glmer 中的缩放预测器无法解决特征值警告；替代优化也没有

Question

lme4 glmer 中的缩放预测器无法解决特征值警告；替代优化也没有

lme4我正在使用 R 中的s函数分析数据（如下）。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)

Run Code Online (Sandbox Code Playgroud)

（注意：我需要将其指定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)

Run Code Online (Sandbox Code Playgroud)

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

“无法评估缩放梯度”
“模型未能收敛：具有 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-unidentabile

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 优化器，使用all_fit()optimx 包中的函数。

这是代码：

# 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)

Run Code Online (Sandbox Code Playgroud)

对于未缩放模型 (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>

Run Code Online (Sandbox Code Playgroud)

对于缩放模型 (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>

Run Code Online (Sandbox Code Playgroud)

数据：

该数据集可通过以下链接获取： https: //www.dropbox.com/s/ud50uatztjq4bh9/20181217%20Surveys%20simplified%20data%20for%20stackX.xlsx ?dl=0

结论及要求：

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

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

Answer 1

Ben*_*ker 6

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

Run Code Online (Sandbox Code Playgroud)

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

归档时间：	7 年前
查看次数：	5742 次
最近记录：	7 年前