Re: [R] bVar slot of lmer objects and standard errors

From: Douglas Bates <dmbates_at_gmail.com>
Date: Sat 31 Dec 2005 - 08:05:20 EST

On 12/29/05, Doran, Harold <HDoran@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@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@stat.math.ethz.ch
> [mailto:r-help-bounces@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@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.@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@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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Dec 31 08:11:05 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:41:44 EST