[R-sig-ME] large dataset - lmer2 vector size specified is too large

[Douglas Bates]

On 3/2/07, Gregor Gorjanc <gregor.gorjanc at bfro.uni-lj.si> wrote:
> Douglas Bates <bates at ...> writes:
> ...
> > but even that formula means that you are estimating 448 fixed effects
> > for which the model matrix is of size 448 * 43680 * 8 bytes (about 150
> > MB).  In addition you are attempting to estimate
>
> Not directly related. Is it really necessary to build model matrix i.e. X.

Necessary? For linear mixed models no - for generalized linear mixed
models yes (think of how the IRLS algorithm works).

Easy to change?  Regretably, no.

> Building X'X directly would be more efficient i.e. only 448 * 448 * 8 bytes.

That wouldn't be sufficient.  As described in the "Implementation"
vignette in the lme4 package, the calculations in lmer2 are based on
the matrix A which is [Z X -y]'[Z X -y].  You can't use just X'X - you
need the crossproducts of all the columns of Z, X and y.

The matrix A is potentially huge but sparse.  Assuming that Florian
really does want to try to estimate a model with (affyID|cephID) then
Z has 224 * 195 = 43680 columns.

In lmer Z'Z is stored as a sparse matrix.  In lmer2 A is stored as a
sparse matrix.

Generating A as a sparse matrix directly from the data or directly
from a model.frame object would result in considerable savings in
memory usage but doing so is far from trivial.  Most useRs are not
aware (nor need they be aware) of the complexity of the operations
involved in taking a linear model formula and a data specification
(and all the arguments like subset, na.action, offset, weights, ...
that are expected to be available in model-fitting functions) and
creating first a model.frame then a model.matrix.  It is one of those
things that "just works" but is incredibly complicated beneath the
surface.

Even the first stage (producing a model frame) is complicated - to the
extent that Thomas Lumley has a report on "The Standard Non-standard
Evaluation in R" dealing with how arguments are evaluated in the
model.frame function.

There are various approaches that could be taken to building A as a
sparse matrix directly from the model specification and a model frame.
 Probably the easiest to implement would be to use the existing
model.matrix code to produce horizontal sections of [Z X -y] and
accumulate A from these sections.  This is still far from trivial
because you either need two passes (one to determine the positions of
the nonzeros and one to evaluate them) or you have to somehow update
the structure of A on the fly.  The other thing that gets tricky is
that the columns of X may change from section to section because not
all levels of a factor may be used in a given section.  It would be
necessary to ensure that the columns from all sections are aligned.