[R-sig-ME] getVarCov for lmer()

[Douglas Bates]

On Feb 1, 2008 11:15 AM, David Afshartous <dafshartous at med.miami.edu> wrote:

> All,

> I've fitted some models using lme() and now I want to check some things by
> fitting them in lmer() (latest version loaded from CRAN).  Under the listed
> methods for lmer() I don't seem to find functionality similar to that of
> getVarCov() to obtain the variance-covariance matrix of the fitted model.
> Does such functionality exist?

I would need to go back and check what getVarCov  produces.  Can you
tell me what you want to be able to calculate?

Be aware that it is not always possible to reproduce all the
capabilities of the nlme package in lme4 because nlme does not handle
(except through egregious kludges) crossed or partially crossed random
effects specifications.  Some of the quantities that can be calculated
for models with a single level of random effects or with multiple
nested levels of random effects cannot even be defined for models with
crossed random effects.

> PS - running R v2.6.0.

> PPS - if it must be computed manually, how does one access information
> stored in slots (e.g., Z', the transpose of the model matrix for the random
> effects)?

If the fitted model is fm then Z' is available as

fm at Zt

It is stored as a compressed, column-oriented sparse matrix of class
"dgCMatrix", defined in the Matrix package.  Linear algebra operations
with such objects should magically do the right thing.  If you are
working with very large models, of course, you will find that you need
to be very careful about exactly what operations are done with Z.  In
many books on mixed models there are formulas for many calculations
based on a generalized least squares representation of the model.
Don't even consider using such formulas except for trivial examples.

I've seen far too many papers where the authors write down an n by n
matrix like V = sigma^2 I + Z Sigma Z' and then proceed to blithely
describe calculations involving the inverse of V.  I usually stop
reading at that point because I know that the author has never tried
to perform the calculation on a non-trivial example.  When n is in the
hundreds of thousands or in the millions you do not even construct n
by n matrices, let alone try to decompose them.  And anyone with
experience in numerical linear algebra knows that you never, ever,
ever expect to evaluate a formula that involves the inverse of a large
matrix.

This is all to say that you could write down an expression involving
solve() applied to some matrix of size n by n derived from Z but you
will be disappointed in the speed.