[R] How to efficiently extract or construct the residual covariance matrix from lme()?

[Andrew Robinson]

Dear R-community,

I'm trying to get the estimated residual covariance matrices from an
lme object.  If we write the model as:

Y = X \beta + Z b + \epsilon

and assume that b ~ N(0, P) and \epsilon ~ N(0, \Sigma), where P is
non-diagonal and \Sigma might have correlation and weights components,
then I'm looking for efficient estimates of

\Sigma

and

ZPZ' + \Sigma

I can find P easily enough, but I'm wondering if there's an easy way
to get at Z and \Sigma.  Also I can move to lmer() if that simplifies
the problem.

Cheers

Andrew
-- 
Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-9763
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
Email: a.robinson at ms.unimelb.edu.au         http://www.ms.unimelb.edu.au