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

Wed Dec 27 22:30:57 CET 2017

This is more likely to get a helpful response if you post on the
r-sig-mixed-models list rather than r-help.

Cheers,
Bert

Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )

On Mon, Dec 25, 2017 at 6:36 PM, Aleksander Główka <aglowka at stanford.edu> wrote:
> 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]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.