New version of lme4 - memory error

Hi Doug,

sorry, it's me having problems again :(

I can install and load the new lme4 package with no trouble, but when
I try to run the examples, I get:
 require(lme4)
Loading required package: lme4
Loading required package: Matrix
Loading required package: lattice
[1] TRUE
sessionInfo()
R version 2.4.1 Patched (2007-01-25 r40572) 
i386-unknown-freebsd6.1 

locale:
C

attached base packages:
[1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
"methods"  
[7] "base"     

other attached packages:
       lme4      Matrix     lattice 
"0.9975-11"  "0.9975-8"   "0.14-16"
fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm1 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy)
Error in as.double(start) : Calloc could not allocate (169499040 of 4)
       memory

Does anyone else find this?  Please let me know what else I can do to
help.

Cheers,

Andrew
Version 0.9975-11 of the lme4 package has been uploaded to CRAN.  The
source package should be available on the mirrors in a day or two and
binary packages should follow soon after.

There are several changes in this release of the package.  The most
important is the availability of a development version of lmer called,
for the time being, lmer2.  At present lmer2 only fits linear mixed
models.  Generalized linear mixed models will be added "soon".
Furthermore there is no mcmcsamp method for a model fit by lmer2.
This deficiency will also be rectified "soon".  Once I have all the
capabilities and methods currently available for lmer also available
for the new representation I will remove the old representation and
rename lmer2 as lmer.

The current version of lmer will continue to be available throughout
the migration process.  You don't have to change anything about your
use of that function unless you want to try the new one.  It would be
a good idea, however, to save the data and the call to lmer in
addition to saving an lmer object, if you so choose, so that you can
recreate the fitted model when the development version becomes the
release version.

The package contains a vignette giving the details of the new implementation.

The reason I am releasing a development version in parallel with the
production version is because I would like feedback from useR's
regarding the development version.  In my experience, testing it
myself and with colleagues whom I visited recently, I have found that
lmer2 is faster and more reliable than the current lmer.  In
particular, on some difficult model fits I have been able to get
substantially better parameter estimates (i.e. the deviance at the
lmer2 estimates is perhaps 4 or 5 lower than that at the lmer
estimates) with lmer2 than I could with lmer.

If you have fit a linear mixed model using lmer and are willing to try
it with lmer2 I would appreciate your telling me if the parameter
estimates are comparable and which fit was faster (use system.time()
to check).  I'm primarily interested in models fit to large data sets
or "difficult" fits.

We have established a new mailing list, R-SIG-mixed-models, for
discussion of R software to fit mixed-effects models, especially lmer.
 See https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models for
information or to subscribe.

I know that I have said this before but this is the last time that I
am going to change the underlying representation.  Really - trust me -
this is the last time.  My theory of software development is expressed
in a line from an old blues song, "you just keep doing it wrong till
you do it right".  I'm convinced that this time I have it right.  That
statement sounds like "famous last words", doesn't it?  :-)

_______________________________________________
R-packages mailing list
R-packages at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-packages

______________________________________________
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.

Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
"Andrew" == Andrew Robinson <A.Robinson at ms.unimelb.edu.au>
    on Sat, 27 Jan 2007 10:27:14 +1100 writes:
Andrew> Hi Doug, sorry, it's me having problems again :(

    Andrew> I can install and load the new lme4 package with no
    Andrew> trouble, but when I try to run the examples, I get:

    >> require(lme4)
    Andrew> Loading required package: lme4
    Andrew> Loading required package: Matrix
    Andrew> Loading required package: lattice
    Andrew> [1] TRUE
    >> sessionInfo()
    Andrew> R version 2.4.1 Patched (2007-01-25 r40572) 
    Andrew> i386-unknown-freebsd6.1 

    Andrew> locale:
    Andrew> C

    Andrew> attached base packages:
    Andrew> [1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
    Andrew> "methods"  
    Andrew> [7] "base"     

    Andrew> other attached packages:
    Andrew> lme4      Matrix     lattice 
    Andrew> "0.9975-11"  "0.9975-8"   "0.14-16" 
    >> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
    >> fm1 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy)
    Andrew> Error in as.double(start) : Calloc could not allocate (169499040 of 4)
    Andrew> memory

Hmm, I can't replicate your problem.  I get

 > system.time( fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
 [1] 0.112 0.001 0.124 0.000 0.000
 > system.time( fm1.2 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy))
 [1] 0.056 0.000 0.056 0.000 0.000
 > summary(fm1.2)
 Linear mixed-effects model fit by REML 
   AIC  BIC logLik MLdeviance REMLdeviance
  1754 1770 -871.8       1752         1744
 Random effects:
  Groups   Name        Variance Std.Dev. Corr  
  Subject  (Intercept) 612.114  24.7409        
	   Days         35.072   5.9222  0.066 
  Residual             654.937  25.5917        
 Number of obs: 180, groups: Subject, 18

 Fixed effects:
	     Estimate Std. Error t value
 (Intercept)  251.405      6.825   36.84
 Days          10.467      1.546    6.77

 Correlation of Fixed Effects:
      (Intr)
 Days -0.138
 > 

    Andrew> Does anyone else find this?  Please let me know what else I can do to
    Andrew> help.

You can type  
        traceback()

after the error (which you should get a habit of doing ;-)
which might be revealing though I doubt it a bit in this case.

Are you sure that you started R as "R --vanilla", i.e.
that you have *not* loaded an .RData of a previous session
accidentally?

Regards,
Martin

    Andrew> Cheers,
    Andrew> Andrew
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)

(I will add it to my list of New Year's Resolutions!)

Andrew
"Andrew" == Andrew Robinson <A.Robinson at ms.unimelb.edu.au>
    on Sat, 27 Jan 2007 10:27:14 +1100 writes:
    Andrew> Hi Doug, sorry, it's me having problems again :(

    Andrew> I can install and load the new lme4 package with no
    Andrew> trouble, but when I try to run the examples, I get:

    >> require(lme4)
    Andrew> Loading required package: lme4
    Andrew> Loading required package: Matrix
    Andrew> Loading required package: lattice
    Andrew> [1] TRUE
    >> sessionInfo()
    Andrew> R version 2.4.1 Patched (2007-01-25 r40572) 
    Andrew> i386-unknown-freebsd6.1 

    Andrew> locale:
    Andrew> C

    Andrew> attached base packages:
    Andrew> [1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
    Andrew> "methods"  
    Andrew> [7] "base"     

    Andrew> other attached packages:
    Andrew> lme4      Matrix     lattice 
    Andrew> "0.9975-11"  "0.9975-8"   "0.14-16" 
    >> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
    >> fm1 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy)
    Andrew> Error in as.double(start) : Calloc could not allocate (169499040 of 4)
    Andrew> memory

Hmm, I can't replicate your problem.  I get

 > system.time( fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
 [1] 0.112 0.001 0.124 0.000 0.000
 > system.time( fm1.2 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy))
 [1] 0.056 0.000 0.056 0.000 0.000
 > summary(fm1.2)
 Linear mixed-effects model fit by REML 
   AIC  BIC logLik MLdeviance REMLdeviance
  1754 1770 -871.8       1752         1744
 Random effects:
  Groups   Name        Variance Std.Dev. Corr  
  Subject  (Intercept) 612.114  24.7409        
	   Days         35.072   5.9222  0.066 
  Residual             654.937  25.5917        
 Number of obs: 180, groups: Subject, 18

 Fixed effects:
	     Estimate Std. Error t value
 (Intercept)  251.405      6.825   36.84
 Days          10.467      1.546    6.77

 Correlation of Fixed Effects:
      (Intr)
 Days -0.138
 > 

    Andrew> Does anyone else find this?  Please let me know what else I can do to
    Andrew> help.

You can type  
        traceback()

after the error (which you should get a habit of doing ;-)
which might be revealing though I doubt it a bit in this case.

Are you sure that you started R as "R --vanilla", i.e.
that you have *not* loaded an .RData of a previous session
accidentally?

Regards,
Martin

    Andrew> Cheers,
    Andrew> Andrew

Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
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.

(I will add it to my list of New Year's Resolutions!)

Andrew

On Sat, Jan 27, 2007 at 04:20:30PM +0100, Martin Maechler wrote:
"Andrew" == Andrew Robinson <A.Robinson at ms.unimelb.edu.au>
    on Sat, 27 Jan 2007 10:27:14 +1100 writes:
    Andrew> Hi Doug, sorry, it's me having problems again :(

    Andrew> I can install and load the new lme4 package with no
    Andrew> trouble, but when I try to run the examples, I get:

    >> require(lme4)
    Andrew> Loading required package: lme4
    Andrew> Loading required package: Matrix
    Andrew> Loading required package: lattice
    Andrew> [1] TRUE
    >> sessionInfo()
    Andrew> R version 2.4.1 Patched (2007-01-25 r40572)
    Andrew> i386-unknown-freebsd6.1

    Andrew> locale:
    Andrew> C

    Andrew> attached base packages:
    Andrew> [1] "stats"     "graphics"  "grDevices" "utils"     "datasets"
    Andrew> "methods"
    Andrew> [7] "base"

    Andrew> other attached packages:
    Andrew> lme4      Matrix     lattice
    Andrew> "0.9975-11"  "0.9975-8"   "0.14-16"
    >> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
    >> fm1 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy)
    Andrew> Error in as.double(start) : Calloc could not allocate (169499040 of 4)
    Andrew> memory

Hmm, I can't replicate your problem.  I get

 > system.time( fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
 [1] 0.112 0.001 0.124 0.000 0.000
 > system.time( fm1.2 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy))
 [1] 0.056 0.000 0.056 0.000 0.000
 > summary(fm1.2)
 Linear mixed-effects model fit by REML
   AIC  BIC logLik MLdeviance REMLdeviance
  1754 1770 -871.8       1752         1744
 Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept) 612.114  24.7409
         Days         35.072   5.9222  0.066
  Residual             654.937  25.5917
 Number of obs: 180, groups: Subject, 18

 Fixed effects:
           Estimate Std. Error t value
 (Intercept)  251.405      6.825   36.84
 Days          10.467      1.546    6.77

 Correlation of Fixed Effects:
      (Intr)
 Days -0.138
 >

    Andrew> Does anyone else find this?  Please let me know what else I can do to
    Andrew> help.

You can type
        traceback()

after the error (which you should get a habit of doing ;-)
which might be revealing though I doubt it a bit in this case.

Are you sure that you started R as "R --vanilla", i.e.
that you have *not* loaded an .RData of a previous session
accidentally?

Regards,
Martin

    Andrew> Cheers,
    Andrew> Andrew
--
Andrew Robinson
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

-------------- next part --------------
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")
[1] 0.51639778 0.09673017 0.00000000
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.

Andrew
Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
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'll respond about how I use gdb in a susequent message.
Andrew
--
Andrew Robinson
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

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 :)

Andrew
Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
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?