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

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: identifying convergence or non-convergence of mixed-effects regression model in lme4 from model output

Bert Gunter-2
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 <[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.