Dear lme4 experts,
Yesterday, I ran the code for two published papers (de Boeck et al.,2011;
de Boeck and Partchev, 2012) on psychometric modeling with glmer in lme4
version 1.1-6 and the vast majority of the models I ran produce convergence
warnings (even the simple ones).
For instance, the basic Rasch model in de Boeck et al. (2011) yields:
## our example data set is bundled with lme4
data("VerbAgg")
## A first example: A Rasch model with fixed item effects and random
+ control=glmerControl(optCtrl=list(maxfun=100000)))
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.308607 (tol = 0.001)
I am a bit puzzeled because, to my knowledge, especially the models for the
VerAgg data (included in lme4) have been checked in many other programs
(also ltm in R) and I heard that glmer produces results that are valid and
consistent with SAS, HLM, ltm, and so on. However, this dataset produces
convergence warnings even though the models are comparatively simple for
psychometric research (basic Rasch and LLTM) and the estimates all seem
reasonable.
I also tried some datasets in the ltm package. Again, convergence warnings
(the mobility data). The estimates are close to what ltm gives me.
Even when I simulate data from a basic Rasch model using the rmvlogis
function in ltm with a couple of extreme item parameters (which occur in
psychometric tests), I also get these warnings. This happens despite glmer
seems to recover the true values very well. The "gradient" errors and the
hessian errors tend to increase with sample size and when I use
binomial("probit") instead of binomial("logit"). It also seems like vector
random effects produce many warnings. optimizer="bobqya" tends to reduce
the warnings but not consistently. To me, it seems like the warnings occur
randomly all over the place. Sometimes glmer "converges" (no warnings) with
one optimizer setting and the values that are very close to the true
values. However, with another optimizer setting, one gets practically
exactly the same estimates and virtually the same likelihood value but a
warning. I really do not understand what is going on.
I had no warnings using version 1.0-5 and version 1.0-6 so this seems to be
a recent problem of lme4?
Is it best to ignore all these convergence warnings for now? Should I
switch back to an older version of lme4 to avoid this problem? Should I
generally avoid using large datasets with lme4?
Many thanks in advance,
Tom
Papers (scripts are online):
De Boeck, P., Bakker, M., Zwitser, R., Nivard, M., Hofman, A., Tuerlinckx,
F., & Partchev, I. (2011). The estimation of item response models with the
lmer function from the lme4 package in R. Journal of Statistical Software,
39, 1-28. http://www.jstatsoft.org/v39/i12
De Boeck, P., & Partchev, I. (2012). IRTrees: Tree-based item response
models of the GLMM family. Journal of Statistical Software, 48, 1-28.
http://www.jstatsoft.org/v48/c01/
Simulated data:
family=binomial,
+ control=glmerControl(optimizer="bobyqa"))
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: value ~ 0 + variable + (1 | person)
Data: simraschlong
AIC BIC logLik deviance df.resid
63887.88 63971.49 -31934.94 63869.88 79991
Random effects:
Groups Name Std.Dev.
person (Intercept) 0.9737
Number of obs: 80000, groups: person, 10000
Fixed Effects:
variableX1 variableX2 variableX3 variableX4 variableX5 variableX6
variableX7 variableX8
3.00598 2.02857 1.01952 0.01191 -1.02146 -1.94484
-2.96534 -9.65211
Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 0.0128742 (tol = 0.001)
[[alternative HTML version deleted]]
Hi,
Sorry for the length of the message below but I also would need some help
and advice.
I am having similar issues (warnings about convergence, may be more) than
other people on this subject but I have been totally lost among all
suggestions and their meanings.
I test the effect of various covriant including inbreeding (O or 1) on
number of occurence of a given behaviour. As it is counting data, I use a
poisson family
The output of my model is :
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 1.26172 (tol = 0.001)
2: In if (resHess$code != 0) { :
the condition has length > 1 and only the first element will be used
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
Analysis of Deviance Table (Type III Wald chisquare tests)
Response: behaviour
Chisq Df Pr(>Chisq)
(Intercept) 1.2814 1 0.257647
inb 2.1574 1 0.141880
inbopp 8.5069 1 0.003538 **
sexratio 3.1415 1 0.076323 .
oppaggre 17.3554 1 3.1e-05 ***
oppmass 0.0165 1 0.897644
duration 358.0516 1 < 2.2e-16 ***
inb:sexratio 7.2522 1 0.007081 **
inbopp:oppmass 8.1315 1 0.004350 **
As usual, I looked at the distribution of residuals: nicely normal.
After that, according to the suggestions, I tested that:
relgrad <- with(test at optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad))
[1] 0.0002681561
-> if this is <0.001 or 0.002, is it enough to assume that the model is ok?
In addition, I tested other optimizer, even if I do not really understand
what it means.
I started with:
testbobyqa <- update(test,control=glmerControl(optimizer="bobyqa"))
[1] 0.001706096
both showed the same warnings and similar results than the first model.
I continued with:
library(optimx)
testnlminb <- update(test,control=glmerControl(optimizer="optimx",
optCtrl=list(method="nlminb")))
testBFGS <- update(test,control=glmerControl(optimizer="optimx",
optCtrl=list(method="L-BFGS-B")))
-> in both case: error
testnlminb:
Error in ans.ret[meth, ] <- c(ans$par, ans$value, ans$fevals, ans$gevals, :
number of items to replace is not a multiple of replacement length
In addition: Warning messages:
1: In optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, :
Parameters or bounds appear to have different scalings.
This can cause poor performance in optimization.
It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
2: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) :
Cholmod warning 'not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 431
3: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) :
Cholmod warning 'not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 431
testBFGS:
Error in eval(expr, envir, enclos) :
(maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate
In addition: Warning messages:
1: In optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, :
Parameters or bounds appear to have different scalings.
This can cause poor performance in optimization.
It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA.
2: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) :
Cholmod warning 'not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 431
3: In pwrssUpdate(pp, resp, tolPwrss, GQmat, compDev, fac, verbose) :
Cholmod warning 'not positive definite' at
file:../Cholesky/t_cholmod_rowfac.c, line 431
4: In optwrap(optimizer, devfun, start, rho$lower, control = control, :
convergence code 9999 from optimx
I then used:
library(nloptr)
## from https://github.com/lme4/lme4/issues/98:
defaultControl <- list(algorithm="NLOPT_LN_BOBYQA",xtol_rel=1e-6,maxeval=1e5)
nloptwrap2 <- function(fn,par,lower,upper,control=list(),...) {
for (n in names(defaultControl))
if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
res <- nloptr(x0=par,eval_f=fn,lb=lower,ub=upper,opts=control,...)
with(res,list(par=solution,
fval=objective,
feval=iterations,
conv=if (status>0) 0 else status,
message=message))
}
test.bobyqa2 <- update(g0.bobyqa,control=glmerControl(optimizer=nloptwrap2))
test.NM2 <- update(g0.bobyqa,control=glmerControl(optimizer=nloptwrap2,
optCtrl=list(algorithm="NLOPT_LN_NELDERMEAD")))
-> similar results than test, testbobyqa, testnelder and warnings for both:
Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 1.88422 (tol = 0.001)
2: In if (resHess$code != 0) { :
the condition has length > 1 and only the first element will be used
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?
relgrad <- with(test.bobyqa2 at optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad))
[1] 1.845328e-05
relgrad <- with(test.NM2 at optinfo$derivs,solve(Hessian,gradient))
max(abs(relgrad))
[1] 1.925339e-05
I also tried that:
getpar <- function(x) c(getME(x,c("theta")),fixef(x))
modList <- list(bobyqa=testbobyqa,NM=testnelder ,
+ bobyqa2=test.bobyqa2,NM2=test.NM2)
ctab <- sapply(modList,getpar)
library(reshape2)
mtab <- melt(ctab)
library(ggplot2)
theme_set(theme_bw())
ggplot(mtab,aes(x=Var2,y=value,colour=Var2))+
+ geom_point()+facet_wrap(~Var1,scale="free")
->I don't really understand wht does this graph mean, if you could explain
to me?
I also looked at this:
likList <- sapply(modList,logLik)
round(log10(max(likList)-likList),1)
bobyqa NM bobyqa2 NM2
-2.3 -2.3 -9.2 -Inf
->What do those values mean?
lbound <- getME(test,"lower")
theta <- getME(test,"theta")
any(lbound==0 & theta<1e-8)
[1] FALSE
-> it means that it is good?
Finally, I tried rescalling my covariables to see if it could be the source
of warnings by using:
datamal$variable.scaled <- datamal$variable/sd(datamal$variable)
and when I introduced the sclaed variables in the model, I received this
message:
Error in FUN(X[[1L]], ...) :
Invalid grouping factor specification, noms:famille
I have tried many things and I am totally lost so if you can help me
clarifying all that, it would be very appreciated.
I hope I wrote down all you needed to understand.
Thanks in advance for your help and sorry again for thlength of the message.
Dear lme4 experts,
Yesterday, I ran the code for two published papers
(de Boeck et al.,2011;
de Boeck and Partchev, 2012) on psychometric modeling with glmer in lme4
version 1.1-6 and the vast majority of the models I
ran produce convergence
warnings (even the simple ones).
[major snippage]
I had no warnings using version 1.0-5 and
version 1.0-6 so this seems to be
a recent problem of lme4?
Is it best to ignore all these convergence warnings for now? Should I
switch back to an older version of lme4 to avoid this problem? Should I
generally avoid using large datasets with lme4?
Many thanks in advance,
Tom
Re-iterating what I may have said earlier:
* these warnings are a result of changes in the WARNING behaviour, not a
change in the fitting algorithm. The only major recent change in the
optimizer has been to move in 1.1-6 from Nelder-Mead to bobyqa for
lmer fits <http://cran.r-project.org/web/packages/lme4/news.html>; this
should isn't relevant to glmer fits and we believe should _improve_
results in almost all cases.
* You should probably ignore the convergence warnings; we now believe
that they are (almost entirely) the result of checking the gradients,
rather than the scaled gradients (i.e. the results of solve(Hessian,grad).
* Version 1.1-7, now on Github (available via devtools::install_github)
and with a macos binary on http://lme4.r-forge.r-project.org/repos ,
tests the scaled gradients and should make the problem go away.
* If you want to double-check, either
* test with lme4.0 (although we
have had reports <http://stackoverflow.com/questions/23662589/
r-reverting-to-lme4-0-and-still-getting-inconsistent-results>
(broken URL!) of 'Downdated XtX' problems -- still chasing that down ...
or * use the code at
https://github.com/lme4/lme4/blob/master/misc/issues/allFit.R
to test with a wide range of optimizers.
Ben Bolker