[R-sig-ME] New version of lme4 - memory error

[Douglas Bates]

On 1/29/07, Andrew Robinson <A.Robinson at ms.unimelb.edu.au> wrote:
> On Sun, Jan 28, 2007 at 12:29:29PM -0600, Douglas Bates wrote:
> > On 1/28/07, Andrew Robinson <A.Robinson at ms.unimelb.edu.au> wrote:
> > >On Sun, Jan 28, 2007 at 11:18:56AM -0600, Douglas Bates wrote:
> > >> On 1/27/07, Andrew Robinson <A.Robinson at ms.unimelb.edu.au> wrote:
> > >> >Thanks, Martin.  It was vanilla.
> > >> >
> > >> >
> > >> >> traceback()
> > >> >4: .Call(mer2_getPars, mer)
> > >> >3: as.double(start)
> > >> >2: nlminb(.Call(mer2_getPars, mer), function(x) .Call(mer2_deviance,
> > >> >       .Call(mer2_setPars, mer, x), as.integer(0)), lower =
> > >ifelse(const,
> > >> >       0, -Inf), control = list(trace = cv$msVerbose, iter.max =
> > >> >       cv$msMaxIter,
> > >> >       rel.tol = abs(0.001/.Call(mer2_deviance, mer, 0))))
> > >> >1: lmer2(Reaction ~ Days + (Days | Subject), sleepstudy)
> > >>
> > >> Thanks for including that traceback, Andrew.  As Martin indicated, it
> > >> is difficult for us to diagnose the problem because we are unable to
> > >> reproduce it.
> > >>
> > >> I can tell you what should be happening at that point in the lmer2
> > >> function.  Perhaps you could use
> > >>
> > >> debug(lmer2)
> > >>
> > >> and check under your system what is happening for you.
> > >>
> > >> As I'm sure you can determine just be looking at this traceback, this
> > >> is the point in lmer2 where the REML or ML criterion is to be
> > >> optimized using the optimizer nlminb.  The object 'mer', of class
> > >> 'mer2', is the internal representation of a mixed-effects model.  It
> > >> is described in the Implementation vignette.
> > >>
> > >> The ST slot is a list whose number of components is equal to the
> > >> number of random effects expressions in the model formula.  In this
> > >> case it should have one component and that component should be a 2x2
> > >> matrix (it's a matrix, not a Matrix from the Matrix package).  The C
> > >> code for mer2_getPars allocates a numeric vector of the correct length
> > >> and fills it out with elements of the lower triangles of these
> > >> matrices starting with the diagonal elements of each matrix.  (The
> > >> diagonal elements are the elements of the scale matrices S and the
> > >> elements of the strict lower triangle are the non-trivial elements of
> > >> the unit triangular matrices T.)
> > >>
> > >> I enclose some of the debugger output from my system.  Once the mer
> > >> object is created I check that the 'dims' slot has the expected
> > >> entries.  In particular the element labeled 'nf' (number of grouping
> > >> factors) should be 1 and the nc slot should be an integer vector of
> > >> length 1 and the first element should be 2.
> > >>
> > >> If the 'dims', 'nc' and 'ST' slots all have the expected contents then
> > >> it would be very mysterious how mer2_getPars managed to fail.  I
> > >> suspect that the problem originates one step back in the call to
> > >> mer2_create.  In either case it may be necessary to go in and examine
> > >> the execution of the compiled code using a symbolic debugger in order
> > >> to exactly what's going on.
> > >
> > >Thanks for your detailed response, Doug.  I ran debug(lmer2) as you
> > >suggested, and as you expected, I got the same results as you did.  Of
> > >course I would like to dig further into this problem.  I have limited
> > >experience using gdb, but I have used it.  If you don't mind providing
> > >a brief set of instructions, I will carry them out.
> >
> > Well, actually, that wasn't what I was expecting.  I thought that
> > something would look wrong in those slots.
> >
> > Can you be more specific about getting the same results as I did?  In
> > particular, did you get the same response to
> >
> > .Call("mer2_getPars", mer, PACKAGE = "lme4")
> >
> > If so, I'm not quite sure how it would fail in the same call
> > immediately afterwards.
>
> I'm sorry, what I told you was wrong.  I get:
>
> Browse[1]>
> debug: mer <- .Call(mer2_create, fl, Zt, t(X), as.double(Y), method ==
>     "REML", nc, cnames, fr$offset, fr$weights)
> Browse[1]>
> debug: const <- unlist(lapply(mer at nc, function(n) rep(1:0, c(n, (n *
>     (n - 1))/2))))
> Browse[1]> mer at ST
> $Subject
>           [,1]       [,2]
> [1,] 0.5163978 0.00000000
> [2,] 0.0000000 0.09673017
>
> Browse[1]> mer at dims
>   nf    n    p    q REML glmm
>    1  180    2   36    1    0
> Browse[1]> mer at nc
> Subject
>       2
> Browse[1]> .Call("mer2_getPars", mer, PACKAGE = "lme4")
> Error: Calloc could not allocate (164788288 of 4) memory

> the last of which is clearly different than you got :)

Congratulations, you get to go in with the symbolic debugger.
Fortunately the C function mer2_getPars is very short and (relatively)
simple.  My guess is that somehow the value of ntot is not being
calculated correctly and the failure is occurring in the call to
allocVector after the first 'for' loop.  For this example the value of
ntot should be 3 at that point.  Other values of importance are nf,
which should be 1, and nc[0] (recall that indexing in C is 0-based)
which should be 2.

Can you check those things for me please?