[R] bVar slot of lmer objects and standard errors

[Douglas Bates]

On 12/29/05, Doran, Harold <HDoran at air.org> wrote:
> Uli:
>
> The graphic in the paper, sometimes called a catepillar plot, must be
> created with some programming as there is (as far as I know) not a
> built-in function for such plots. As for the contents of bVar you say
> the dimensions are 2,2,28 and there are two random effects and 28
> schools. So, from what I know about your model, the third dimension
> represents the posterior covariance matrix for each of your 28 schools
> as Spencer notes.
>
> For example, consider the following model
> > library(Matrix)
> > library(mlmRev)
> > fm1 <- lmer(math ~ 1 + (year|schoolid), egsingle)
>
> Then, get the posterior means (modes for a GLMM)
> > fm1 at bVar$schoolid
>
> These data have 60 schools, so you will see ,,1 through ,,60 and the
> elements of each matrix are posterior variances on the diagonals and
> covariances in the off-diags (upper triang) corresponding to the
> empirical Bayes estimates for each of the 60 schools.
>
> , , 1
>
>            [,1]         [,2]
> [1,] 0.01007129 -0.001272618
> [2,] 0.00000000  0.004588049

I'd have to go back and check but I think that these are the upper
triangles of the symmetric matrix (as Spencer suggested) that are the
conditional variance-covariance matrices of the two-dimensional 
random effects for each school up to a scale factor.  That is, I think
each face needs to be multiplied by s^2 to get the actual
variance-covariance matrix.

>
>
> Does this help?
>
> Harold
>
>
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Spencer Graves
> Sent: Thursday, December 29, 2005 6:58 AM
> To: Ulrich Keller
> Cc: r-help
> Subject: Re: [R] bVar slot of lmer objects and standard errors
>
>           Have you received a satisfactory reply to this post?  I
> haven't seen one.  Unfortunately, I can't give a definitive answer, but
> I can offer an intelligent guess.  With luck, this might encourage
> someone who knows more than I do to reply.  If not, I hope these
> comments help you clarify the issue further, e.g., by reading the source
> or other references.
>
>           I'm not not sure, but I believe that
> lmertest1 at bVar$schoolid[,,i] is the upper triangular part of the
> covariance matrix of the random effects for the i-th level of schoolid.
>   The lower triangle appears as 0, though the code (I believe) iterprets
> it as equal to the upper triangle.  More precisely, I suspect it is
> created from something that is stored in a more compact form, i.e.,
> keeping only a single copy of the off-diagonal elements of symmetric
> matrices.  I don't seem to have access to your "nlmframe", so I can't
> comment further on those specifics.  You might be able to clarify this
> by reading the source code.  I've been sitting on this reply for several
> days without finding time to do more with it, so I think I should just
> offer what I suspect.
>
>           The specifics of your question suggest to me that you want to
> produce something similar to Figure 1.12 in Pinheiro and Bates (2000)
> Mixed-Effects Models in S and S-Plus (Springer).  That was produced from
> an "lmList" object, not an "lme" object, so we won't expect to get their
> exact answers.  Instead, we would hope to get tighter answers available
> from pooling information using "lme";  the function "lmList" consideres
> each subject separately with no pooling.  With luck, the answers should
> be close.
>
>           I started by making a local copy of the data:
>
> library(nlme)
> OrthoFem <- Orthodont[Orthodont$Sex=="Female",]
>
>           Next, I believe to switch to "lme4", we need to quit R
> completely and restart.  I did that.  Then with the following sequence
> of commands I produced something that looked roughly similar to the
> confidence intervals produced with Figure 1.12:
>
> library(lme4)
> fm1OrthF. <- lmer(distance~age+(age|Subject), data=OrthoFem)
>
> fm1.s <-  coef(fm1OrthF.)$Subject
> fm1.s.var <- fm1OrthF. at bVar$Subject
> fm1.s0.s <- sqrt(fm1.s.var[1,1,])
> fm1.s0.a <- sqrt(fm1.s.var[2,2,])
> fm1.s[,1]+outer(fm1.s0.s, c(-2,0,2))
> fm1.s[,2]+outer(fm1.s0.a, c(-2,0,2))
>
>           hope this helps.
>           Viel Glueck.
>           spencer graves
>
> Ulrich Keller wrote:
>
> > Hello,
> >
> > I am looking for a way to obtain standard errors for emprirical Bayes
> estimates of a model fitted with lmer (like the ones plotted on page 14
> of the document available at
> http://www.eric.ed.gov/ERICDocs/data/ericdocs2/content_storage_01/000000
> 0b/80/2b/b3/94.pdf).
>
>
> Harold Doran mentioned
> (http://tolstoy.newcastle.edu.au/~rking/R/help/05/08/10638.html)
> that  the posterior modes' variances can be found in the bVar slot of
> lmer objects. However, when I fit e.g. this model:
> >
> > lmertest1<-lmer(mathtot~1+(m_escs_c|schoolid),hlmframe)
> >
> > then lmertest1 at bVar$schoolid is a three-dimensional array with
> dimensions (2,2,28).
> The factor schoolid has 28 levels, and there are random effects for the
> intercept and m_escs_c, but what does the third dimension correspond to?
> In other words, what are the contents of bVar, and how can I use them to
> get standard errors?
> >
> > Thanks in advance for your answers and Merry Christmas,
> >
> > Uli Keller