[R-sig-ME] Wrong convergence warnings with glmer in lme4 version 1.1-6??

UG uriel.gelin at usherbrooke.ca
Wed May 14 16:13:43 CEST 2014


Tom Davis <tomd792 at ...> writes:

> 
> 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
> person effects
> > m1 = glmer(r2 ~ 0 + item + (1|id), data=VerbAgg, family="binomial",
> + 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:
> 
> > library("reshape")
> > library("ltm")
> > library("lme4")
> > set.seed("12345")
> >
> > simrasch<-data.frame(rmvlogis(200, cbind(seq(-5, 5, 0.5), 1)))
> > rasch(simrasch,IRT.param=F)
> 
> Call:
> rasch(data = simrasch, IRT.param = F)
> 
> Coefficients:
>     X1      X2      X3      X4      X5      X6      X7      X8      X9
> X10     X11     X12
> 17.342   5.157   5.157   3.441   3.141   2.100   2.193   1.969   1.484
> 0.717   0.054  -0.347
>    X13     X14     X15     X16     X17     X18     X19     X20
> X21       z
> -0.874  -1.415  -2.103  -2.696  -3.065  -3.152  -4.453  -4.216  -5.175
> 1.091
> 
> Log.Lik: -1329.186
> 
> >
> > simrasch$person = rownames(simrasch)
> > simraschlong = melt(simrasch, id = "person")
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> family=binomial)
> 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
>  2702.484  2842.026 -1329.242  2658.484      4178
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 1.091
> Number of obs: 4200, groups:  person, 200
> Fixed Effects:
>  variableX1   variableX2   variableX3   variableX4   variableX5
> variableX6   variableX7
>    19.04000      5.14239      5.13514      3.43274      3.13736
> 2.10085      2.19335
>  variableX8   variableX9  variableX10  variableX11  variableX12
> variableX13  variableX14
>     1.97086      1.48698      0.71828      0.05206     -0.35255
> -0.88099     -1.42318
> variableX15  variableX16  variableX17  variableX18  variableX19
> variableX20  variableX21
>    -2.11070     -2.70069     -3.06758     -3.15545     -4.44743
> -4.21083     -5.16615
> Warning messages:
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge with max|grad| = 0.135068 (tol = 0.001)
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> 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
>  2702.482  2842.025 -1329.241  2658.482      4178
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 1.09
> Number of obs: 4200, groups:  person, 200
> Fixed Effects:
>  variableX1   variableX2   variableX3   variableX4   variableX5
> variableX6   variableX7
>    17.75458      5.14416      5.14417      3.43797      3.13991
> 2.10444      2.19676
>  variableX8   variableX9  variableX10  variableX11  variableX12
> variableX13  variableX14
>     1.97412      1.48907      0.72049      0.05389     -0.34882
> -0.87841     -1.42030
> variableX15  variableX16  variableX17  variableX18  variableX19
> variableX20  variableX21
>    -2.10751     -2.69767     -3.06471     -3.15126     -4.44429
> -4.20819     -5.16356
> 
> > simrasch<-data.frame(rmvlogis(1000, cbind(seq(-5, 5, 0.5), 1)))
> > rasch(simrasch,IRT.param=F)
> 
> Call:
> rasch(data = simrasch, IRT.param = F)
> 
> Coefficients:
>     X1      X2      X3      X4      X5      X6      X7      X8      X9
> X10     X11     X12
>  5.195   4.422   3.868   3.599   3.192   2.558   2.211   1.605   0.984
> 0.559   0.003  -0.398
>    X13     X14     X15     X16     X17     X18     X19     X20
> X21       z
> -0.985  -1.590  -1.994  -2.401  -2.876  -3.411  -3.815  -4.683  -4.832
> 1.014
> 
> Log.Lik: -6874.631
> 
> >
> > simrasch$person = rownames(simrasch)
> > simraschlong = melt(simrasch, id = "person")
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> family=binomial)
> 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
> 13793.707 13968.657 -6874.853 13749.707     20978
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 1.014
> Number of obs: 21000, groups:  person, 1000
> Fixed Effects:
>  variableX1   variableX2   variableX3   variableX4   variableX5
> variableX6   variableX7
>    5.182206     4.415991     3.871713     3.599747     3.191841
> 2.561370     2.219830
>  variableX8   variableX9  variableX10  variableX11  variableX12
> variableX13  variableX14
>    1.611983     0.992610     0.564159     0.003716    -0.398537
> -0.987771    -1.591732
> variableX15  variableX16  variableX17  variableX18  variableX19
> variableX20  variableX21
>   -1.998437    -2.400455    -2.870716    -3.408691    -3.806151
> -4.672642    -4.813643
> Warning message:
> In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge with max|grad| = 0.907091 (tol = 0.001)
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> 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
> 13793.696 13968.646 -6874.848 13749.696     20978
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 1.014
> Number of obs: 21000, groups:  person, 1000
> Fixed Effects:
>  variableX1   variableX2   variableX3   variableX4   variableX5
> variableX6   variableX7
>    5.184774     4.414100     3.862784     3.594478     3.190096
> 2.559954     2.214893
>  variableX8   variableX9  variableX10  variableX11  variableX12
> variableX13  variableX14
>    1.610161     0.988363     0.562450     0.002823    -0.400700
> -0.990108    -1.595035
> variableX15  variableX16  variableX17  variableX18  variableX19
> variableX20  variableX21
>   -1.997929    -2.403493    -2.875794    -3.408200    -3.809483
> -4.674454    -4.822606
> Warning message:
> In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge with max|grad| = 0.00124419 (tol = 0.001)
> >
> >
> >
> 
> > set.seed("12345")
> > simrasch<-data.frame(rmvlogis(200, cbind(c(seq(-3, 3, 1),10), 1)))
> > rasch(simrasch,IRT.param=F)
> 
> Call:
> rasch(data = simrasch, IRT.param = F)
> 
> Coefficients:
>      X1       X2       X3       X4       X5       X6       X7
> X8        z
>   3.226    1.820    1.461    0.044   -0.706   -1.494   -2.771  -10.454
> 0.801
> 
> Log.Lik: -647.099
> 
> >
> > simrasch$person = rownames(simrasch)
> > simraschlong = melt(simrasch, id = "person")
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> family=binomial)
> 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
> 1312.7956 1361.1955 -647.3978 1294.7956      1591
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 0.7881
> Number of obs: 1600, groups:  person, 200
> Fixed Effects:
> variableX1  variableX2  variableX3  variableX4  variableX5  variableX6
> variableX7  variableX8
>    3.21214     1.82080     1.46473     0.04478    -0.71031    -1.49811
> -2.76154   -27.10659
> Warning message:
> In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> 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
> 1312.7956 1361.1955 -647.3978 1294.7956      1591
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 0.7881
> Number of obs: 1600, groups:  person, 200
> Fixed Effects:
> variableX1  variableX2  variableX3  variableX4  variableX5  variableX6
> variableX7  variableX8
>    3.21213     1.82080     1.46473     0.04478    -0.71032    -1.49811
> -2.76154   -19.56371
> 
> > set.seed("12345")
> > simrasch<-data.frame(rmvlogis(10000, cbind(c(seq(-3, 3, 1),10), 1)))
> > rasch(simrasch,IRT.param=F)
> 
> Call:
> rasch(data = simrasch, IRT.param = F)
> 
> Coefficients:
>     X1      X2      X3      X4      X5      X6      X7      X8       z
>  3.018   2.026   1.011   0.012  -1.013  -1.941  -2.977  -9.696   0.986
> 
> Log.Lik: -31914.71
> 
> >
> > simrasch$person = rownames(simrasch)
> > simraschlong = melt(simrasch, id = "person")
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> family=binomial)
> 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
>  63888.47  63972.08 -31935.24  63870.47     79991
> Random effects:
>  Groups Name        Std.Dev.
>  person (Intercept) 0.9739
> Number of obs: 80000, groups:  person, 10000
> Fixed Effects:
> variableX1  variableX2  variableX3  variableX4  variableX5  variableX6
> variableX7  variableX8
>    3.01884     2.03562     1.02945     0.02299    -1.01066    -1.93548
> -2.95234    -9.44503
> Warning messages:
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge with max|grad| = 25.6865 (tol = 0.001)
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
> > glmer(value ~ 0 + variable + (1|person), data=simraschlong,
> 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?

> summary(test)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) ['glmerMod']
 Family: poisson ( log )
Formula: behaviour~ inbreeding+ opponent inbreeding + sexratio + opponent
aggressiveness+ opponent mass + duration of the test +  
    inbreeding:sexratio + opponent inbreeding:opponent mass+ (1 |
family/IDfocal)
   Data: datamal

     AIC      BIC   logLik deviance df.resid 
  1469.3   1507.5   -723.6   1447.3      228 

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.211 -1.223 -0.279  0.847  4.491 

Random effects:
 Groups       Name        Variance Std.Dev.
 ID:family (Intercept) 0.25751  0.5075  
 family      (Intercept) 0.04598  0.2144  
Number of obs: 239, groups: ID:family, 63; family, 19

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.4300096  0.3798771   1.132  0.25765    
inb             -0.6387831  0.4348942  -1.469  0.14188    
inbopp           1.2304585  0.4218720   2.917  0.00354 ** 
sexratio        -0.7260132  0.4096143  -1.772  0.07632 .  
oppaggre         0.0017798  0.0004272   4.166  3.1e-05 ***
oppmass          0.0010126  0.0078721   0.129  0.89764    
duration         0.0065125  0.0003442  18.922  < 2e-16 ***
inb:sexratio     1.7354628  0.6444371   2.693  0.00708 ** 
inbopp:oppmass  -0.0499209  0.0175064  -2.852  0.00435 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
               (Intr) inb    inbopp sexrat oppaggr oppmass dura  inb:sx
inb            -0.598                                                 
inbopp         -0.241 -0.004                                          
sexratio       -0.761  0.686 -0.007                                   
oppaggre       -0.121  0.021 -0.116  0.013                            
oppmass        -0.545  0.048  0.482  0.052 -0.101                     
duration       -0.079 -0.053 -0.013 -0.019 -0.331 -0.095              
inb:sexrati     0.499 -0.900  0.004 -0.665 -0.039 -0.048  0.100       
inbopp:oppmass  0.198  0.012 -0.987  0.012  0.154 -0.436 -0.001 -0.012

> Anova(test, type="III")
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"))
> max(abs(relgrad))
[1] 0.01506629
testnelder <- update(test,control=glmerControl(optimizer="Nelder_Mead"))
> max(abs(relgrad))
[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.



More information about the R-sig-mixed-models mailing list