title: Scaling predictors in lme4 glmer doesn't resolve eigenvalue warnings; neither does alternative optimization
body: <p>I am analysing data (included below) using <code>lme4</code>'s <code>glmer</code> function in R. 
  The model I am building consists of a Poisson-distributed response variable (<code>obs</code>), one random factor (<code>area</code>), one continuous offset (<code>duration</code>), five continuous fixed effects (<code>can_perc</code>, <code>can_n</code>, <code>time</code>, <code>temp</code>, <code>cloud_cover</code>) and one binomial fixed effect factor (<code>burnt</code>).
  Before fitting the model I checked for collinearity and removed any collinear variables.</p>
  
  <p>The initial model is:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>(Note: I need to specify the <code>na.action</code> as 'na.fail' as I want to <code>dredge()</code> the model later and this is required for that.)</p>
  
  <p>Running the model gives the following warning:</p>
  
  <blockquote>
    <p>"Hessian is numerically singular: parameters are not uniquely determined"</p>
  </blockquote>
  
  <p>In similar variations of the model, I have also received the warning:</p>
  
  <blockquote>
    <p>"In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
    Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?"</p>
  </blockquote>
  
  <p>From my limited understanding of the advice here <a href="https://rdrr.io/cran/lme4/man/troubleshooting.html" rel="nofollow noreferrer">https://rdrr.io/cran/lme4/man/troubleshooting.html</a> and elsewhere, both of these warnings reflect a similar issue, of the Hessian (inverse curvature matrix) having a large eigenvalue, indicating that (within numerical tolerances) the surface is completely flat in some direction.
  Based on the advice in the warnings and link, I rescaled all of the continuous predictor variables using <code>scale()</code>. I also scaled the offset variable (I tried both with and without scaling this one). The model with scaled predictor variables is here:</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>However I have not yet escaped the eigenvalues! The scaled model gives two warnings:</p>
  
  <blockquote>
    <p>"unable to evaluate scaled gradient" <br>
    "Model failed to converge: degenerate  Hessian with 1 negative eigenvalues"</p>
  </blockquote>
  
  <p>I have searched a lot online and couldn't find another case/solution to how to deal with eigenvalue problems once the predictors have been scaled, other than checking that the model hasn't been misspecified.</p>
  
  <h1>Attempts to address warnings/improve optimization:</h1>
  
  <p>Based on these pages/documents: 
  <a href="https://cran.r-project.org/web/packages/lme4/lme4.pdf" rel="nofollow noreferrer">https://cran.r-project.org/web/packages/lme4/lme4.pdf</a></p>
  
  <p><a href="https://stats.stackexchange.com/questions/164457/r-glmer-warnings-model-fails-to-converge-model-is-nearly-unidentifiable">https://stats.stackexchange.com/questions/164457/r-glmer-warnings-model-fails-to-converge-model-is-nearly-unidentifiable</a> </p>
  
  <p><a href="https://rdrr.io/cran/lme4/man/isSingular.html" rel="nofollow noreferrer">https://rdrr.io/cran/lme4/man/isSingular.html</a></p>
  
  <p><a href="https://rdrr.io/cran/lme4/man/troubleshooting.html" rel="nofollow noreferrer">https://stats.stackexchange.com/questions/242109/model-failed-to-converge-warning-in-lmer</a></p>
  
  <p>and others,</p>
  
  <p>I have:</p>
  
  <ol>
  <li><p>checked the model specifications and data for mistakes (none that I can see - have I missed something?) </p></li>
  <li><p>checked for singularity with <code>is_singular(x, tol = 1e-05)</code> (somehow this function call evolved from <code>isSingular()</code> to current form?): all models give FALSE.</p></li>
  <li><p>checked convergence measure with <code>converge_ok(q2, tolerance = 0.001)</code>: All models give FALSE, unless I substantially increase the tolerance; however they do vary considerably in their convergence measure.</p></li>
  <li><p>tried different optimizers/model estimation methods as follows:  </p>
  
  <ul>
  <li>a) <code>glmerControl(optimizer = "bobyqa") and glmerControl(optimizer ="Nelder_Mead")</code>  </li>
  <li>b) <code>glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb'))</code></li>
  <li>c) bobyqa, Nelder_Mead, optimx.nlminb, optimx.L-BFGS-B, nloptwrap.NLOPT_LN_NELDERMEAD, nloptwrap.NLOPT_LN_BOBYQA, and nmkbw optimisers, using the <code>all_fit()</code> function from optimx package.  </li>
  </ul></li>
  </ol>
  
  <p>Here is the code:</p>
  
  <pre><code># 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)
  </code></pre>
  
  <h1>Output from the above code, for the unscaled model (q1):</h1>
  
  <pre><code>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
  &lt;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&gt;
  
  $`optimx.L-BFGS-B`
  &lt;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&gt;
  
  $nloptwrap.NLOPT_LN_NELDERMEAD
  &lt;simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite&gt;
  
  $nloptwrap.NLOPT_LN_BOBYQA
  &lt;simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite&gt;
  
  $nmkbw.
  &lt;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&gt;
  </code></pre>
  
  <h1>Output from the above code, for the scaled model (q2):</h1>
  
  <pre><code>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
  &lt;simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite&gt;
  
  $nloptwrap.NLOPT_LN_BOBYQA
  &lt;simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite&gt;
  
  $nmkbw.
  &lt;simpleError in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev,     grpFac = fac, maxit = maxit, verbose = verbose): Downdated VtV is not positive definite&gt;
  </code></pre>
  
  <h1>Data:</h1>
  
  <p>The dataset is available at this link:
  <a href="https://www.dropbox.com/s/ud50uatztjq4bh9/20181217%20Surveys%20simplified%20data%20for%20stackX.xlsx?dl=0" rel="nofollow noreferrer">https://www.dropbox.com/s/ud50uatztjq4bh9/20181217%20Surveys%20simplified%20data%20for%20stackX.xlsx?dl=0</a></p>
  
  <h1>Conclusion and request:</h1>
  
  <p>It looks to me that none of these alternative methods of optimisation have succeeded either; in fact some of then seem to have raised other warnings/errors which would take me down another rabbit hole.</p>
  
  <p>Can anyone advise how I could progress with the fitting of these models? 
  It is not my intent for these to be the final models, but rather to dredge them and then select optimal/top models from the different alternative subset models.</p>