[R] identifying convergence or non-convergence of mixed-effects regression model in lme4 from model output

[Aleksander Główka]

Hi R community!

I've fitted three mixed-effects regression models to a thousand 
bootstrap samples (case-resampling regression) using the lme4 package in 
a custom-built for-loop. The only output I saved were the inferential 
statistics for my fixed and random effects. I did not save any output 
related to the performance to the machine learning algorithm used to fit 
the models (REML=FALSE). After running all the simulations, I got about 
two dozen messages of this kind:

25: In checkConv(attr(opt, "derivs"), opt$par, ctrl = 
control$checkConv,  ... :
   Model failed to converge with max|grad| = 4.49732 (tol = 0.002, 
component 1)
26: In checkConv(attr(opt, "derivs"), opt$par, ctrl = 
control$checkConv,  ... :
   Model is nearly unidentifiable: large eigenvalue ratio
  - Rescale variables?

Since I only get the error messages after all the computations have been 
executed, looking at these error messages is not helpful, because, as 
you can see, they don't allow me to identify which bootstrap sample the 
non-converged model was fit to they are referring to. Since I don't 
think I can recover this information from the already run simulations, I 
need to modify my information extractor code to include convergence 
information and rerun the simulations.

My question is the following: which attribute in the summary of a 
mixed-effects model in lme4 allows me to check if the model has 
converged or not? What value would the parameter corresponding to that 
attribute have to have in order for me to conclude the model has not 
converged?

Here are my current extractor functions for fixed and random effects.

lmer.data.extract = function(lmer.mod, name=deparse(substitute(lmer.mod))){

   #extract predictor names & create data frame, attach other cols to 
new data frame
   mod.data = data.frame(summary(lmer.mod)$coefficients[,1])
   names(mod.data) = "estimate"
   mod.data$std.error = as.numeric(summary(lmer.mod)$coefficients[,2]) 
#std errors
   mod.data$df = as.numeric(summary(lmer.mod)$coefficients[,3]) #degrees 
of freedom
   mod.data$t.val = as.numeric(summary(lmer.mod)$coefficients[,4]) #t-values
   mod.data$p.val = as.numeric(summary(lmer.mod)$coefficients[,5]) 
#p-values

   #extract AIC, BIC, logLik, deviance df.resid
   mod.data$AIC = as.numeric(summary(lmer.mod)$AIC[1])
   mod.data$BIC = as.numeric(summary(lmer.mod)$AICtab[2][1])
   mod.data$logLik = as.numeric(summary(lmer.mod)$AICtab[3][1])
   mod.data$deviance = as.numeric(summary(lmer.mod)$AICtab[4][1])
   mod.data$df.resid = as.numeric(summary(lmer.mod)$AICtab[5][1])

   #add number of datapoints
   mod.data$N = as.numeric(summary(incr.best.m)$devcomp$dims[1])

   #add model name
   mod.data$model = name

   return(mod.data)

}

lmer.ranef.data.extract = function(lmer.mod, 
name=deparse(substitute(lmer.mod))){

   #extract random effect variance, standard error, correlations between 
slope and intercept
   mod.data.ef = as.data.frame(VarCorr(lmer.mod))

   mod.data.ef$n.subj = as.numeric(summary(lmer.mod)$ngrps[1]) #number 
of subjects
   mod.data.ef$n.item = as.numeric(summary(lmer.mod)$ngrps[2]) #number 
of items

   #add number of datapoints
   mod.data.ef$N = as.numeric(summary(incr.best.m)$devcomp$dims[1])

   #add model name
   mod.data.ef$model = name

   return(mod.data.ef)

}

I'm also including the structure of an example model that did converge 
(but I can I tell from the output?).

List of 18
  $ methTitle   : chr "Linear mixed model fit by maximum likelihood  
\nt-tests use  Satterthwaite approximations to degrees of freedom"
  $ objClass    : atomic [1:1] lmerMod
   ..- attr(*, "package")= chr "lme4"
  $ devcomp     :List of 2
   ..$ cmp : Named num [1:10] 176.85 59.09 95.43 3.84 99.27 ...
   .. ..- attr(*, "names")= chr [1:10] "ldL2" "ldRX2" "wrss" "ussq" ...
   ..$ dims: Named int [1:12] 1742 1742 10 1732 94 4 1 2 0 0 ...
   .. ..- attr(*, "names")= chr [1:12] "N" "n" "p" "nmp" ...
  $ isLmer      : logi TRUE
  $ useScale    : logi TRUE
  $ logLik      :Class 'logLik' : -65 (df=15)
  $ family      : NULL
  $ link        : NULL
  $ ngrps       : Named num [1:2] 36 29
   ..- attr(*, "names")= chr [1:2] "subj" "item"
  $ coefficients: num [1:10, 1:5] 7.00546 0.04234 -0.00258 0.09094 
-0.00804 ...
   ..- attr(*, "dimnames")=List of 2
   .. ..$ : chr [1:10] "(Intercept)" "FreqABCD.log.std" 
"LogitABCD.neg.log.std" "MIABCD.neg.log.std" ...
   .. ..$ : chr [1:5] "Estimate" "Std. Error" "df" "t value" ...
  $ sigma       : num 0.239
  $ vcov        :Formal class 'dpoMatrix' [package "Matrix"] with 5 slots
   .. ..@ x       : num [1:100] 1.15e-03 3.40e-05 -5.12e-05 2.18e-05 
3.65e-06 ...
   .. ..@ Dim     : int [1:2] 10 10
   .. ..@ Dimnames:List of 2
   .. .. ..$ : chr [1:10] "(Intercept)" "FreqABCD.log.std" 
"LogitABCD.neg.log.std" "MIABCD.neg.log.std" ...
   .. .. ..$ : chr [1:10] "(Intercept)" "FreqABCD.log.std" 
"LogitABCD.neg.log.std" "MIABCD.neg.log.std" ...
   .. ..@ uplo    : chr "U"
   .. ..@ factors :List of 1
   .. .. ..$ correlation:Formal class 'corMatrix' [package "Matrix"] 
with 6 slots
   .. .. .. .. ..@ sd      : num [1:10] 0.0339 0.0519 0.013 0.0439 
0.0068 ...
   .. .. .. .. ..@ x       : num [1:100] 1 0.0194 -0.1162 0.0147 0.0158 ...
   .. .. .. .. ..@ Dim     : int [1:2] 10 10
   .. .. .. .. ..@ Dimnames:List of 2
   .. .. .. .. .. ..$ : chr [1:10] "(Intercept)" "FreqABCD.log.std" 
"LogitABCD.neg.log.std" "MIABCD.neg.log.std" ...
   .. .. .. .. .. ..$ : chr [1:10] "(Intercept)" "FreqABCD.log.std" 
"LogitABCD.neg.log.std" "MIABCD.neg.log.std" ...
   .. .. .. .. ..@ uplo    : chr "U"
   .. .. .. .. ..@ factors :List of 1
   .. .. .. .. .. ..$ Cholesky:Formal class 'Cholesky' [package 
"Matrix"] with 5 slots
   .. .. .. .. .. .. .. ..@ x       : num [1:100] 1 0 0 0 0 0 0 0 0 0 ...
   .. .. .. .. .. .. .. ..@ Dim     : int [1:2] 10 10
   .. .. .. .. .. .. .. ..@ Dimnames:List of 2
   .. .. .. .. .. .. .. .. ..$ : NULL
   .. .. .. .. .. .. .. .. ..$ : NULL
   .. .. .. .. .. .. .. ..@ uplo    : chr "U"
   .. .. .. .. .. .. .. ..@ diag    : chr "N"
  $ varcor      :List of 2
   ..$ subj: num [1, 1] 0.0273
   .. ..- attr(*, "dimnames")=List of 2
   .. .. ..$ : chr "(Intercept)"
   .. .. ..$ : chr "(Intercept)"
   .. ..- attr(*, "stddev")= Named num 0.165
   .. .. ..- attr(*, "names")= chr "(Intercept)"
   .. ..- attr(*, "correlation")= num [1, 1] 1
   .. .. ..- attr(*, "dimnames")=List of 2
   .. .. .. ..$ : chr "(Intercept)"
   .. .. .. ..$ : chr "(Intercept)"
   ..$ item: num [1:2, 1:2] 0.00417 0.000484 0.000484 0.00289
   .. ..- attr(*, "dimnames")=List of 2
   .. .. ..$ : chr [1:2] "(Intercept)" "FreqABCD.log.std"
   .. .. ..$ : chr [1:2] "(Intercept)" "FreqABCD.log.std"
   .. ..- attr(*, "stddev")= Named num [1:2] 0.0646 0.0538
   .. .. ..- attr(*, "names")= chr [1:2] "(Intercept)" "FreqABCD.log.std"
   .. ..- attr(*, "correlation")= num [1:2, 1:2] 1 0.139 0.139 1
   .. .. ..- attr(*, "dimnames")=List of 2
   .. .. .. ..$ : chr [1:2] "(Intercept)" "FreqABCD.log.std"
   .. .. .. ..$ : chr [1:2] "(Intercept)" "FreqABCD.log.std"
   ..- attr(*, "sc")= num 0.239
   ..- attr(*, "useSc")= logi TRUE
   ..- attr(*, "class")= chr "VarCorr.merMod"
  $ AICtab      : Named num [1:5] 159.7 241.6 -64.8 129.7 1727
   ..- attr(*, "names")= chr [1:5] "AIC" "BIC" "logLik" "deviance" ...
  $ call        : language lme4::lmer(formula = RT.log ~ 
FreqABCD.log.std + LogitABCD.neg.log.std + MIABCD.neg.log.std + 
AS.data$freq.sub.PC1 +      AS.data$freq.sub.PC2 + AS.data$freq.sub.PC3 
+ AS.data$freq.sub.PC4 + block + nletter.std + (1 | subj) +  ...
  $ residuals   : Named num [1:1742] 0.713 0.498 -0.361 -0.101 2.594 ...
   ..- attr(*, "names")= chr [1:1742] "1" "2" "3" "4" ...
  $ fitMsgs     : chr(0)
  $ optinfo     :List of 7
   ..$ optimizer: chr "bobyqa"
   ..$ control  :List of 1
   .. ..$ iprint: int 0
   ..$ derivs   :List of 2
   .. ..$ gradient: num [1:4] 9.81e-06 -5.34e-06 -1.60e-05 7.06e-05
   .. ..$ Hessian : num [1:4, 1:4] 245.9 28.5 3.3 -13.7 28.5 ...
   ..$ conv     :List of 2
   .. ..$ opt : int 0
   .. ..$ lme4: list()
   ..$ feval    : int 107
   ..$ warnings : list()
   ..$ val      : num [1:4] 0.6919 0.2705 0.0314 0.223
  - attr(*, "class")= chr "summary.merMod"

I'd appreciate any advice you may have!

Thank you,

Aleksander Główka
PhD Candidate
Department of Linguistics
Stanford University
**

	[[alternative HTML version deleted]]