extracting coefficients from lmer2 model

I meant to cc: the list on this reply

---------- Forwarded message ----------
From: Douglas Bates <bates at stat.wisc.edu>
Date: Feb 1, 2007 1:44 PM
Subject: Re: [R-sig-ME] extracting coefficients from lmer2 model
To: Colin Beale <c.beale at macaulay.ac.uk>

Hi all,

Some of you may have seen my recent message to the r-help list
concerning a modification to lmer code that caused R to crash. Happily,
the same modifications to the lmer2 code do not cause this crash so as
far as I am concerned this is a major improvement to the original. Now,
however, I'm not sure how to get at the bits of the fitted object that I
want to play with to post-process the results as the methods aren't yet
implemented for this structure - I need, for example, the coefficients
of the fixed and random effects and the full variance-covariance matrix
of fixed _and_ random effects (together!). I'm sure there must be ways
around this issue, since I can see that there are slots with names that
look very suitable at least to hold the fixed and random effects - but I
can't work out how to get them. And I'm not sure if the full
variance-covariance matrix is going to be in there anyway (looks like it
might be for lmer objects).

I recommend the vignette in the lme4 package that describes the lmer2
representation.

In answer to some of your specific questions:

fixef is straightforward for the mer2 representation

The ranef slot is currently populated when you call for fixef but the
values in that slot are the random effects with the "spherical" prior
density (independent, constant-variance).  These are named b* in the
vignette.   I need to add code to undo the permutation and to undo the
orthogonalization transformation defined by S and T.  It is not
incredibly difficult but I will need some quiet time to consider
exactly what sequence of operations needs to be performed and then to
write and debug the code.

The variance-covariance matrix for the random effects and the fixed
effects can be derived from the inverse of L but do you really want to
create that object?  It has the potential of being very large and very
slow to calculate.  Earlier today there was discussion of the need to
set control = list(niterEM = 0, gradient = FALSE) in a call to lmer
that fits a complicated model to a large data set.  This reason that
the ECME iterations and the gradient calculations are so slow (on the
order of hundreds or thousands of evaluations of the deviance) is
because those calculations require part of the inverse of L.  It's
even simpler than that, all you need to do is to calculate some groups
of columns in the inverse and store some values derived from them.
You don't ever need to store the inverse itself - just these values
that are derived from groups of columns in the inverse.

Calculating the inverse of L won't be a problem for simple models and
small data sets but it will be a problem for multiple, non-nested
grouping factors with many levels.  Furthermore the sparse matrix
methods take care to minimize the amount of storage in L.  That's what
the fill-reducing permutation is designed to do.  However, the inverse
of L can be a mess.

You should consider carefully exactly what you need from that
variance-covariance matrix.  For example, if you are going to evaluate
a quadratic form using this matrix then don't calculate the matrix.
Write the quadratic form using the factorization and solve the system
involving L or L-transpose then collect an inner product.

For example:

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm2 <- lmer2(Reaction ~ Days + (Days|Subject), sleepstudy)

using:

 slotNames(fm1)
 slotNames(fm2)

shows me that there are slots with the ranef and fixef in both classes
of object - I can use fixef() on both objects to get the fixed effects,
but ranef() only works on fm1 at the moment - how else can I get these
objects directly? And there are obviously many more slots in the fm1
object than the fm2 object - can anyone point out if any of these slots
are the full variance-covariance matrix I need, containing both fixed
and random effects and if so, how I can get at it (or if not, how I
might go about creating it?).

Thanks in advance,

Colin


Dr. Colin Beale
Spatial Ecologist
The Macaulay Institute
Craigiebuckler
Aberdeen
AB15 8QH
UK

Tel: 01224 498245 ext. 2427
Fax: 01224 311556
Email: c.beale at macaulay.ac.uk


--
Please note that the views expressed in this e-mail are those of the
sender and do not necessarily represent the views of the Macaulay
Institute. This email and any attachments are confidential and are
intended solely for the use of the recipient(s) to whom they are
addressed. If you are not the intended recipient, you should not read,
copy, disclose or rely on any information contained in this e-mail, and
we would ask you to contact the sender immediately and delete the email
from your system. Thank you.
Macaulay Institute and Associated Companies, Macaulay Drive,
Craigiebuckler, Aberdeen, AB15 8QH.


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models