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.

In particular, the problems occur in models for paths in the analysis that link two variables measured only at plot-level, so the data fed into the model have plot-level variables repeated for each individual observation within the same plot, and random effect of plot to account for the grouping of trees within plots.  Outside of this path analysis context, I would condense these data to the plot-level values and run the model with the unit of replication as the plot, without the random effects.

A simplified subset of the data (about 1/3 the individuals), with each row an individual tree -- can be found here: https://dl.dropboxusercontent.com/u/9541528/path.data.example.R


    plot species survival      size     env1       env2      env3
71 BS195    LIDE        0 -1.258929 1.536898 -0.8162132 0.2175846
72 BS195    LIDE        0 -1.258929 1.536898 -0.8162132 0.2175846
73 BS196    QUAG        1  1.888022 1.261792  0.2997178 0.4836638
74 BS196    UMCA        1  1.182543 1.261792  0.2997178 0.4836638
75 BS196    UMCA        0  1.173050 1.261792  0.2997178 0.4836638

Survival is binary 1/0; size has been measured at the tree level and scaled/centered. Env1, env2, env3 are measured at the plot level and have been scaled/centered.

Among others, I get convergence failures because of very large max grad values, degenerate Hessian with 1 (or more) negative eigenvalues, something about "bobyqa -- a trust region step failed to reduce q", or unidentifiable models with very large eigenvalues. 

An example model giving me problems:
lmer(env1 ~ env2 + env3 + (1|plot) + (1|species), data=data)
Linear mixed model fit by REML ['merModLmerTest']
Formula: env1 ~ env2 + env3 + (1 | plot) + (1 | species)
   Data: 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)
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. 

           - - - - - - - -

lme() from the nlme package does not give me the same errors, but ends up with very different estimates..  I have not been using lme() because some of my models require binomial errors for a logistic regression approach, and because I have limited understanding of how to cross two random effects in lme() syntax.  Setting that aside for a moment, here are outputs from two models, as an example:

> summary(lme(env1 ~ env2 + env3, random=~1|plot, data=data))
Linear mixed-effects model fit by REML
 Data: data 
        AIC       BIC   logLik
  -22563.92 -22544.59 11286.96

Random effects:
 Formula: ~1 | plot
        (Intercept)     Residual
StdDev:    1.239896 3.326088e-16

Fixed effects: env1 ~ env2 + env3 
                 Value  Std.Error  DF   t-value p-value
(Intercept)  0.5808555 0.05079791 332 11.434633  0.0000
env2        -0.0679534 0.04492343  21 -1.512650  0.1453
env3         0.1970801 0.06227939  21  3.164452  0.0047
 Correlation: 
     (Intr) env2 
env2 0.136       
env3 0.457  0.019

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.3379242 -0.6675848  0.0000000  0.8344810  2.0027545 

Number of Observations: 356
Number of Groups: 24 


Compared to lmer() from lme4:


> 	lmer(env1 ~ env2 + env3 + (1|plot), data=data)
Linear mixed model fit by REML ['merModLmerTest']
Formula: env1 ~ env2 + env3 + (1 | plot)
   Data: data
REML criterion at convergence: -6469.918
Random effects:
 Groups   Name        Std.Dev. 
 plot     (Intercept) 2.434e-01
 Residual             1.276e-05
Number of obs: 356, groups: plot, 24
Fixed Effects:
(Intercept)         env2         env3  
     0.1078       0.2882      -0.2411  
Warning messages:
1: In optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp),  :
  convergence code 3 from bobyqa: bobyqa -- a trust region step failed to reduce q
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0239058 (tol = 0.002)

           - - - - - - - -

If I set my model equal to g0 and run the code that Ben suggested in response to another posting earlier today (to scale the estimated gradient by the curvature/standard errors), I get NAs because of the negative eigenvalue in the hessian (as far as I understand):

> gg <- g0@optinfo$derivs$grad
> hh <- g0@optinfo$derivs$Hessian
> vv <- sqrt(diag(solve(hh/2)))
Warning message:
In sqrt(diag(solve(hh/2))) : NaNs produced
> summary(abs(gg*vv))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  1.184   1.184   1.184   1.184   1.184   1.184       1 

I would appreciate any suggestions on how to understand what is happening, what checks I should run to assess whether these are false positives, what the implications of a singular model fit are, etc... 

Thank you so very much for any guidance you might provide.
Best
Margaret Metz


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Margaret R. Metz
Assistant Project Scientist
Department of Plant Pathology
One Shields Ave.
University of California
Davis, CA 95616
	[[alternative HTML version deleted]]

