[R-sig-ME] lme4 convergence warnings and confirmatory path analysis with hierarchical data

[Ben Bolker]

Margaret Metz <mrmetz at ...> writes:

> 
> Hello all -- 
 
> I am attempting a confirmatory path analysis that uses lme4 for
> linear mixed effects models with hierarchically structured data
> (following methods in Shipley 2009 Ecology --
> http://www.esajournals.org/doi/pdf/10.1890/08-1034.1). The data are
> about tree survival with predictors measured at the level of the
> individual tree or at the level of the plot within which the trees
> are found.
 
> My models are triggering errors and warnings that have been the
> subject of quite a few recent posts on this list.  I have tried some
> of the checks that Ben Bolker has generously suggested, but I do not
> know how to interpret the outcomes.  The errors/warnings have also
> been changing as new versions of lme4 are available.  I would
> appreciate help to know whether it is possible for me to resolve
> these errors or whether I should seed a different analysis approach.
 

  Short answer -- I haven't gotten very far into this yet.
It did inspire me to post a utility function for trying out
all known optimizers:

https://raw.githubusercontent.com/lme4/lme4/master/misc/issues/allFit.R

When I try this on your data,

## setwd("~/R/pkgs/lme4git/lme4/misc/issues")
source("path.data.example.R")
library(lme4)
fm1 <- lmer(env1 ~ env2 + env3 + (1|plot) + (1|species), data=data)
source("allFit.R") ## URL above
rr <- allFit(fm1)
rr2 <- rr[!sapply(rr,inherits,"try-error")]
sapply(rr2,fixef)  ## bobyqa, nlminb, L-BFGS-B all give similar fixed effects
sapply(rr2,function(x) unlist(VarCorr(x)))   ## and RE variances
lapply(rr2,function(x) x at optinfo$conv$lme4$messages)  ## all complain
lapply(rr2,function(x) x at optinfo$warnings)

But my results are quite different from yours -- are we using the
same subset of the data?

> REML criterion at convergence: -9807.396
> Random effects:
>  Groups   Name        Std.Dev. 
>  plot     (Intercept) 2.555e-01
>  species  (Intercept) 4.179e-12
>  Residual             8.358e-08
> Number of obs: 356, groups: plot, 24; species, 9
> Fixed Effects:
> (Intercept)         env2         env3  
>     0.02395      0.56535     -0.16874  
> Warning messages:
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge with max|grad| = 3795.17 (tol = 0.002)

  This message goes away in my local branch, where I scale the
gradient differently

> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
 
> Some models run with the full dataset give errors that
> "Downdated VtV is not positive definite."  Some model
> fits are singular, I think, but I do not fully understand
> the implications of that for proceeding.

  You know what it means, right? (Zero variances/perfect correlations
in the estimated random effects variance-covariance matrix).  I haven't
followed Shipley 2009, but if everything _else_ about your data
were reliable, and if your focus was on the fixed-effect estimates,
I might not worry about it too much.  Or you could remove those
terms from the model.  It really depends a bit more on the context
of your analysis.