[R] bVar slot of lmer objects and standard errors

Fri Dec 30 23:14:44 CET 2005

Hi, Doug:

	  Thanks.  Perhaps the easiest way to check would be to find an example 
with answers you believe to be correct and compare what this gives you 
with the presumed "correct" answers.  I don't have such an example 
handy, but Ulrich Keller, who originated this thread, might.  If not, 
maybe he can find another way to check the answer, e.g., via Monte Carlo.

	  Best Wishes,
	  Spencer Graves

Douglas Bates wrote:

> 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
> 
> 
> ______________________________________________
> R-help at 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

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915