R help - Dec 2005 - bVar slot of lmer objects and standard errors

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/0000000b/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

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/0000000b/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

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


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

______________________________________________
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