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