[R-sig-ME] Structure of an mer object and steps in the computation within lmer

[Douglas Bates]

This is in response to David Duffy's query about creating specialized
Z matrices for lmer.  I am responding in a separate message to
generate a new thread, rather than having the discussion hidden in an
misleading topic.

All of lmer, glmer and nlmer (which now works, by the way) now create
the same class of object, class "mer".  The declaration of the class,
in the file lme4/R/AllClass.R, has some comments regarding each of the
slots.  They are arranged so that the original data, model matrices,
etc. come first, then the slots that vary during the optimization.

X, the model matrix for the fixed effects, is a (dense) matrix in the
usual orientation (rows corresponding to observations).  Zt is the
transpose of the model matrix for the random effects stored as a
sparse matrix of class "dgCMatrix" defined in the Matrix package. The
y slot is, unsurprisingly, the response.  The prior weights and the
offset, if used, are stored as numeric (double precision, not integer)
slots of length n.  If they are not used they must still be numeric
vectors but they are of length 0.  This is for the S4 class
definitions.  In S3 these would typically be NULL if not used but in
S4 they must always be the same class so numeric(0) means "not used".

Random effects are associated with different terms in the model
formula.  The terms not only define sections of the model matrix Z but
also determine the structure of the variance-covariance of the random
effects.  The Gp slot is the "groups pointer" that indicates which
elements of the random effects vector are associated with each random
effects term.  To get the association of levels and covariates to
random effects within the terms, think of the result of
ranef(fittedModel).  It is a list of matrices, one for each term.  The
number of rows of the ith matrix is the number of levels of the
grouping factor for the ith term.

Try fitting models with different numbers of terms and different
expressions on the left hand side of the | to see how this all plays
out.  For a fitted model fm the expression

all.equal(unname(unlist(ranef(fm))), fm at ranef)

should be TRUE.  The form of Gp is an integer vector whose first
element is 0 and successive elements are the cumulative sum of the
length of the contents of the ranef matrix for each term.  If that
explanation is too obscure, try a few examples to see what I mean.

The elements of the dims slot and deviance slots are documented at the
beginning of the lmer.c page at lme4.r-forge.r-project.org/doxygen/

The ST slot is a list (one element per term) of square matrices.  The
ith element is of size ncol(ranef(fm)[[i]]).  They contain the
elements of the relative variance-covariance matrices of the random
effects coded in a particular way.  You don't have to know the details
if you just .Call the C function "mer_ST_initialize" to create initial
values.

The slot A is another sparse matrix and, except for nonlinear models,
it has the same strucuture as Zt so you can initialize it to that
value.  There is a C function "mer_create_L" to create the L slot from
the Zt matrix.  The L slot is the sparse Cholesky factor of a matrix
derived from Z'Z.  It is a CHMfactor object from the Matrix package.

There is another sparse matrix C that is used for generalized linear
mixed models and for nonlinear mixed models but in slightly different
forms.  For a GLMM it has the same structure as A so only the
numerical values are stored in the Cx slot.  For a NLMM the entire
sparse matrix is stored in the Cm slot.

The simple form of the optimization is to create an mer object, run
the validity checks on it

validObject(ans)

then call mer_finalize on the object with an integer value of verbose
(0 for no output during iterations, 1 for a 1-line summary every
iteration, -1 for even more output from GLMMs and NLMMs).

The details of exactly what happens inside mer_finalize are different
for linear mixed models than for other types of mixed models.  For an
LMM the fixed effects and the common scale parameter (i.e. sigma) are
profiled out of the optimization so only the parameters in ST are
varied.  The ST parameters are updated with "mer_ST_setPars" and the
profiled deviance for both ML and REML estimation is evaluated by
"mer_update_RX".

For GLMMs or NLMMs you must set both the ST parameters and the fixef
slot, then call "mer_update_u" to get the conditional modes of the
transformed random effects u.  For a GLMM the deviance is evaluated
from the deviance residuals in "mer_update_dev".

I hope this helps.  I'm sure there will be questions, assuming that
anyone is crazy enough and brave enough to want to mess with this
stuff.