Aleksander Główka
2017-Dec-26 02:36 UTC
[R] identifying convergence or non-convergence of mixed-effects regression model in lme4 from model output
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]]
Bert Gunter
2017-Dec-27 21:30 UTC
[R] identifying convergence or non-convergence of mixed-effects regression model in lme4 from model output
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.