[R-sig-ME] lmer nonconvergent: care to run and explain?

Douglas Bates bates at stat.wisc.edu
Fri Oct 16 19:44:09 CEST 2015 

Instead of saying  "the unscaled one still fails", Paul, it would be better
to say that the unscaled model produces an warning about convergence.  In
some cases it is easy to verify that the algorithm has converged; in other
cases it isn't.

As an aside, when we were developing nlme I often needed to alert people to
the other case - declaring convergence when in fact the algorithm had not
converged.  Many people told us how superior the nonmem software was to the
nlme package because nonmem always converged, which was true in the sense
that it always declared convergence but it did that whether or not it was
anywhere close to the optimum.

If you look at the trace of the parameters for your unscaled case, the
fiinal values are effectively on the boundary.  Those two parameters are
ratios of standard deviations and must be non-negative.  The first value is
the standard deviation of the random effects for the intercept divided by
the standard deviation of the per-observation random noise term.  There two
standard deviations are on the same scale, which is the scale of the
response.  A value of 5e-7 is essentially zero.  In the MixedModels package
for Julia I check after convergence has been declared whether there are
small values in the parameter vector and, if so, set them to zero and check
that the objective has not been increased.  Almost all the time setting
these small values to zero decreases the objective so, instead of reporting
this ratio as 5e-7 it reports the more easily recognized value of zero.

Regarding the warning I don't think the calculations used here are
meaningful for convergence on the boundary.  There are different conditions
to be met when you converge to an interior point and when you converge on
the boundary.

Colin Longhurst and i have been checking the performance of various
optimizers used in R and in Julia for fitting linear mixed-effects models.
The approach is simple - we use whatever examples we can find in the
literature and fit them using different optimizers that are available in
the two languages.  The results are in a package called "Timings" available
at

https://github.com/Stat990-033/Timings

The results are in a directory called JSON of the installed package.  (Not
surprisingly, the results are saved in Javascript Object Notation,
abbreviated JSON).

I would like to add your example and its timings if i may.  To be sure that
we are fitting the same model I have

> str(dat)'data.frame':	15419 obs. of  8 variables:
 $ Mind: int  1 1 1 1 1 1 1 1 1 1 ...
 $ Iind: num  1 2 3 4 5 6 7 8 9 10 ...
 $ y1  : num  -37.037 -24.231 -0.254 -41.583 26.973 ...
 $ x1  : num  75.3 120 103.2 106.5 89.6 ...
 $ x2  : num  199 202 214 232 226 ...
 $ x3  : num  128 156 190 158 103 ...
 $ y2  : num  -38.45 -25.64 -1.67 -43 25.56 ...

$ y3 : num -31.74 -14.96 7.52 -33.51 33.54 ...

and the result of fitting

lmer(formula = y3 ~ 1 + x1 + x2 + x3 + (1 | Mind) + (0 + x1 | Mind),
data = dat, verbose = 3L)


ends with

318:     143253.09: 1.22600e-06 0.00975153
319:     143253.09:  0.00000 0.00975152
320:     143253.09: 1.12044e-07 0.00975136
321:     143253.09: 1.82473e-07 0.00975160
322:     143253.09:  0.00000 0.00975153
At return
eval: 322 fn:      143253.09 par:  0.00000 0.00975153Warning
message:In checkConv(attr(opt, "derivs"), opt$par, ctrl =
control$checkConv,  :
  Model is nearly unidentifiable: very large eigenvalue
 - Rescale variables?


(Notice that in my case the first value does come out as 0.  Also the
number of iterations are different.)


I am using

other attached packages:

[1] lme4_1.1-10  Matrix_1.2-2 MASS_7.3-44


loaded via a namespace (and not attached):

[1] minqa_1.2.4


In summary, I don't think there is anything to worry about from those
warnings.  You can go back and refit the model omitting the (1|Mind)
term and you should get the same value of the REML criterion.


On Fri, Oct 16, 2015 at 4:57 AM Thierry Onkelinx <thierry.onkelinx at inbo.be>
wrote:

> Dear Paul,
>
> I can reproduce the problem with lme4_1.1-9 and Matrix_1.2-3.
>
> Lowering the SE of the noise solved the covergence issues. So I expect that
> a low signal to noise ratio is causing the false convergence.
>
> Best regards,
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be no more
> than asking him to perform a post-mortem examination: he may be able to say
> what the experiment died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer does not
> ensure that a reasonable answer can be extracted from a given body of data.
> ~ John Tukey
>
> 2015-10-15 16:49 GMT+02:00 Paul Johnson <pauljohn32 at gmail.com>:
>
> > On Wed, Oct 14, 2015 at 7:19 PM, David Duffy
> > <David.Duffy at qimrberghofer.edu.au> wrote:
> > > Is it because Mind is not declared a factor?
> > >
> > > dat$Mind <- factor(dat$Mind)
> > > lmer( y3s ~ x1s + x2s + x3s + (x1s | Mind), data = dat)
> > > Linear mixed model fit by REML ['lmerMod']
> > > Formula: y3s ~ x1s + x2s + x3s + (x1s | Mind)
> > >    Data: dat
> > > REML criterion at convergence: 33271.78
> > > Random effects:
> > >  Groups   Name        Std.Dev. Corr
> > >  Mind     (Intercept) 0.6837
> > >           x1s         0.1373   1.00
> > >  Residual             0.6997
> > > Number of obs: 15419, groups:  Mind, 100
> > > Fixed Effects:
> > > (Intercept)          x1s          x2s          x3s
> > >   -0.003039     0.117080    -0.163680     0.055911
> > >
> > >
> > > | David Duffy (MBBS PhD)
> > > | email: David.Duffy at qimrberghofer.edu.au  ph: INT+61+7+3362-0217 fax:
> > -0101
> > > | Genetic Epidemiology, QIMR Berghofer Institute of Medical Research
> > > | 300 Herston Rd, Brisbane, Queensland 4006, Australia  GPG 4D0B994A
> >
> > Interesting. I think there's a ghost in the machine.
> >
> > I see no benefit.
> >
> > The unscaled one still fails:
> >
> > > dat$Mind <- as.factor(dat$Mind)
> > > mm3 <- lmer( y3 ~ x1 + x2 + x3 + (1|Mind) + (0 + x1 | Mind), data=dat,
> > verbose=3)
> >
> > ...
> > 350:     143253.09: 9.88456e-06 0.00975148
> > 351:     143253.09: 9.88260e-06 0.00975168
> > 352:     143253.09: 9.69049e-06 0.00975152
> > 353:     143253.09: 9.44279e-06 0.00975152
> > 354:     143253.09: 9.61561e-06 0.00975152
> > At return
> > eval: 354 fn:      143253.09 par: 9.44279e-06 0.00975152
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
> >   Model is nearly unidentifiable: very large eigenvalue
> >  - Rescale variables?
> >
> > And the scaled one is still a failure on my system
> >
> > dat[ , c("x1s", "x2s", "x3s")] <- lapply(dat[, c("x1", "x2", "x3")],
> scale)
> > dat$y3s <- scale(dat$y3)
> > mm3 <- lmer( y3s ~ x1s + x2s + x3s + (1 | Mind) + (0 + x1 | Mind),
> >             data=dat, verbose = 3)
> >
> > rho:  2.0e-07 eval: 186 fn:      33271.8 par: 0.00000 0.00975124
> > 187:     33271.783:  0.00000 0.00975151
> > 188:     33271.783: 5.40934e-07 0.00975153
> > 189:     33271.783: 1.08231e-06 0.00975153
> > 190:     33271.783: 6.32312e-07 0.00975153
> > At return
> > eval: 190 fn:      33271.783 par: 5.40934e-07 0.00975153
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
> >   Model is nearly unidentifiable: very large eigenvalue
> >  - Rescale variables?
> >
> >
> > Since you don't get errors, I wonder if it means you are using
> > different version of R, or lmer?  Can we see your sessionInfo for
> > comparison:
> >
> >
> > > sessionInfo()
> > R version 3.2.2 (2015-08-14)
> > Platform: x86_64-pc-linux-gnu (64-bit)
> > Running under: Ubuntu 15.04
> >
> > locale:
> >  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
> >  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
> >  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
> >  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
> >  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> >
> > attached base packages:
> > [1] stats     graphics  grDevices utils     datasets  methods   base
> >
> > other attached packages:
> > [1] lme4_1.1-10  Matrix_1.2-2 MASS_7.3-44
> >
> > loaded via a namespace (and not attached):
> > [1] minqa_1.2.4     compiler_3.2.2  tools_3.2.2     Rcpp_0.12.1
> > [5] splines_3.2.2   nlme_3.1-122    grid_3.2.2      nloptr_1.0.4
> > [9] lattice_0.20-33
> >
> >
> > --
> > Paul E. Johnson
> > Professor, Political Science        Director
> > 1541 Lilac Lane, Room 504      Center for Research Methods
> > University of Kansas                 University of Kansas
> > http://pj.freefaculty.org              http://crmda.ku.edu
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]

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