[R] lme() doesn't converge on IGF example

[Douglas Bates]

On 6/13/07, David Daniel <ddaniel at nmsu.edu> wrote:
> Running the Chapter 4 examples in Pinheiro & Bates' "Mixed-Effects
> Models in S and S-PLUS" (2000), I get a message that the default
> optimizer doesn't converge, but using "optim" for the optimizer
> results in convergence:
>
> > > library(nlme)
> > > fm1IGF.lis <- lmList(IGF)
> > > fm1IGF.lme <- lme(fm1IGF.lis)
> > Error in lme.formula(fixed = conc ~ age, data = IGF, random = list
> > (Lot = c(-0.741604809797216,  :
> >       nlminb problem, convergence error code = 1; message = iteration
> > limit reached without convergence (9)
> > >
> > > fm1IGF.lme <- lme(fm1IGF.lis, control= list(opt="optim"))
>
> I wouldn't have expected the default optimizer to not work with an
> example from this text.  Not knowing anything about the optimizers,
> I'm wondering if this is expected or known behavior, or if there are
> tips for getting it to converge other than changing optimizers?

That model fit corresponds to a singular variance-covariance matrix
for the random effects (notice that the correlation is -1).  The way
that the model was written in lme this corresponds to an infinite
value of one of the parameters so it is actually an advantage that the
nlminb optimizer doesn't declare convergence.

In the lmer2 function from the lme4 package the model is defined in
such a way that the singular variance-covariance matrix corresponds to
a value of zero for one of the parameters that is constrained to be
nonnegative.   Try

library(lme4)
data(IGF, package = "nlme")
fm1IGF.lmer <- lmer2(conc ~ age + (age|Lot), IGF, control =
list(msVerbose = TRUE))

and you will see that the second parameter is exactly zero at
convergence (or, at least it is on my amd_64 Linux system).

Currently I do not flag this for the user in the "show" method for the
fitted model.  I should.  Depending on your point of view such a
fitted model is either a boundary case or not a legitimate mixed
model.

The ability to converge to a singular model is actually the big
difference between the lmer and the lmer2 functions in the lme4
package.  For the lmer2 function the model is expressed in such a way
that the log-likelihood or the REML criterion can be evaluated for
singular variance-covariance matrices.  Furthermore such evaluations
approach the boundary evaluations smoothly.  The lmer function
evaluates the log-likelihood using the precision matrix (i.e. the
inverse of the variance-covariance) which, by definition, cannot be
evaluated when the variance-covariance matrix is singular.

>
> nlme Version:       3.1-80
>
> > > R.Version()
> > $platform
> > [1] "i386-apple-darwin8.9.1"
> >
> > $arch
> > [1] "i386"
> >
> > $os
> > [1] "darwin8.9.1"
> >
> > $system
> > [1] "i386, darwin8.9.1"
> >
> > $status
> > [1] ""
> >
> > $major
> > [1] "2"
> >
> > $minor
> > [1] "5.0"
> >
> > $year
> > [1] "2007"
> >
> > $month
> > [1] "04"
> >
> > $day
> > [1] "23"
> >
> > $`svn rev`
> > [1] "41293"
> >
> > $language
> > [1] "R"
> >
> > $version.string
> > [1] "R version 2.5.0 (2007-04-23)"
>
> ----------------------------------
> David Daniel
> Associate Professor
> University Statistics Center
> New Mexico State University
>
> ddaniel at nmsu.edu
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> 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.
>