bVar slot of lmer objects and standard errors

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

bVar slot of lmer objects and standard errors

Ulrich Keller
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@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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Spencer Graves
          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/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@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
>
> ______________________________________________
> [hidden email] 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

[hidden email]
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Doran, Harold
In reply to this post by Ulrich Keller
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


Does this help?

Harold


-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] 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
>
> ______________________________________________
> [hidden email] 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

[hidden email]
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Douglas Bates
On 12/29/05, Doran, Harold <[hidden email]> 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: [hidden email]
> [mailto:[hidden email]] 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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Spencer Graves
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 <[hidden email]> 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: [hidden email]
>>[mailto:[hidden email]] 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
>
>
> ______________________________________________
> [hidden email] 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

[hidden email]
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Ulrich Keller
In reply to this post by Douglas Bates
Hello,

I'm sorry to resurrect this thread that I started almost two months ago.
I've been pretty busy since I posted my question and the issue is not
that high on my priority list. Thanks to all those who replied, and I
hope I can tickle your interest again.

As a reminder, my question was how one can extract the conditional
posterior variance of a random effect from the bVar slot of an lmer
model. Thanks to your answers, I now understand that I have to use the
diagonal elements of the conditional matrices. However, I am not quite
sure what this means:

Douglas Bates wrote:
> 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.

What is s^2? Where can I find it in the lmer object? I tried reading the
source, but gave up fairly quickly. Thanks in advance for your replies,
and this time I promise I'll be more responsive.


Uli

My original post:

> 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@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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Spencer Graves
          I shall provide herewith an example of what I believe is "the
conditional posterior variance of a random effect".  I hope someone more
knowledgeable will check this and provide a correction if it's not
correct.  I say this in part because my simple rationality check (below)
didn't match as closely as I had hoped.

          I don't have easy access to the "hlmframe" of your example, so I used
the example in the "lmer" documentation:

library(lme4)
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
Formula: Reaction ~ Days + (Days | Subject)
    Data: sleepstudy
       AIC      BIC    logLik MLdeviance REMLdeviance
  1753.628 1769.593 -871.8141   1751.986     1743.628
Random effects:
  Groups   Name        Variance Std.Dev. Corr
  Subject  (Intercept) 612.090  24.7405
           Days         35.072   5.9221  0.066
  Residual             654.941  25.5918
# of obs: 180, groups: Subject, 18
<snip>

          I think we want the "Residual Variance" of 654.941 (Std.Dev. of
25.5918).  To get this, let's try "VarCorr":

(vC.fm1 <- VarCorr(fm1))
$Subject
2 x 2 Matrix of class "dpoMatrix"
             (Intercept)     Days
(Intercept)   612.09032  9.60428
Days            9.60428 35.07165

attr(,"sc")
[1] 25.59182

          Our desired 25.5918 is 'attr(,"sc")' here.  We get that as follows:

(s. <- attr(vC.fm1, "sc"))
[1] 25.59182

          Next, by studying 'str(fm1)' and earlier emails you cite, we get the
desired conditional posterior covariance matrix as follows:

 > (condPostCov <- (s.^2)*fm1@bVar$Subject)
, , 1
          [,1]       [,2]
[1,] 145.7051 -21.444432
[2,]   0.0000   5.312275
<snip>
, , 18
          [,1]       [,2]
[1,] 145.7051 -21.444432
[2,]   0.0000   5.312275

          This actually gives us only the upper triangular portion of the
desired covariance matrices.  The following will make them all symmetric:

condPostCov[2,1,] <- condPostCov[1,2,]

          As a check, I computed the covariance matrix of the estimated random
effects. To get this, I reviewed str(ranef(fm1)), which led me to the
following:

  var(ranef.fm1@.Data[[1]])
             (Intercept)     Days
(Intercept)   466.38503 31.04874
Days           31.04874 29.75939

          These numbers are all much larger than "condPostCov".  However, I
believe this must be a random bounce -- unless it's a deficiency in my
understanding (in which case, I hope someone will provide a correction).

          Ulrich:  Would you mind checking this either with a published example
or a Monte Carlo and reporting the results to us?

          Viel Glück,
          spencer graves
         
Ulrich Keller wrote:

> Hello,
>
> I'm sorry to resurrect this thread that I started almost two months ago.
> I've been pretty busy since I posted my question and the issue is not
> that high on my priority list. Thanks to all those who replied, and I
> hope I can tickle your interest again.
>
> As a reminder, my question was how one can extract the conditional
> posterior variance of a random effect from the bVar slot of an lmer
> model. Thanks to your answers, I now understand that I have to use the
> diagonal elements of the conditional matrices. However, I am not quite
> sure what this means:
>
> Douglas Bates wrote:
>
>>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.
>
>
> What is s^2? Where can I find it in the lmer object? I tried reading the
> source, but gave up fairly quickly. Thanks in advance for your replies,
> and this time I promise I'll be more responsive.
>
>
> Uli
>
> My original post:
>
>>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@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
>
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: bVar slot of lmer objects and standard errors

Ulrich Keller
Vielen Dank, Spencer.

I could not find a published example where both the original data and
conditional posterior variances were available. Instead, I toyed around
a little with artificial data, and the (pretty pathetic) result below is
the closest I came to "Monte Carlo-ing". I'm afraid I lack the
statistical and R skills to do it properly (in a reasonable amount of
time). Still, the result looks about right.

I started with a some simple artificial data:
 > y<-rep(100,100)+rep(c(-10,-5,2,15),c(20,30,40,10))+rnorm(100,0,15)
 > f<-factor(rep(c(1,2,3,4),c(20,30,40,10)))
 > fy.df<-data.frame(f,y)

To which I then fitted a simple model:
 > fy.lmer<-lmer(y~1 + (1 | f), data=fy.df)

Now to get independent estimates of the conditional posterior variances,
I used the bootstrap, re-fitting the model to 100 resampled dataframes
(sampled independently for each group) and storing the resulting random
effect estimates in a matrix:
 > boots<-matrix(NA,100,4)
 > colnames(boots)<-levels(fy.df$f)
 > for (i in 1:nrow(boots)) {
+   resampled.df<-do.call("rbind", lapply(split(fy.df, fy.df$f),#
+     function(x) x[sample(nrow(x), nrow(x), replace=TRUE),]))
+   boots[i,]<-ranef(update(fy.lmer,data=resampled.df))[1]$f
+ }

Finally, I compared the bootstrap estimates to those provided by lmer:
 > s<-attr(VarCorr(fy.lmer),"sc")
 > rbind(fy.lmer@bVar$f[1,,]*s^2,apply(boots,2,var))
            1        2        3        4
[1,] 7.908634 5.448064 4.155261 14.42237
[2,] 5.935841 4.501380 5.367947 14.04685

Which looks quite good. I tried this for several artificial datasets,
and the result always looked ok.



Spencer Graves wrote:

>       I shall provide herewith an example of what I believe is "the
> conditional posterior variance of a random effect".  I hope someone
> more knowledgeable will check this and provide a correction if it's
> not correct.  I say this in part because my simple rationality check
> (below) didn't match as closely as I had hoped.
>
>       I don't have easy access to the "hlmframe" of your example, so I
> used the example in the "lmer" documentation:
>
> library(lme4)
> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
> Formula: Reaction ~ Days + (Days | Subject)
>    Data: sleepstudy
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  1753.628 1769.593 -871.8141   1751.986     1743.628
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr
>  Subject  (Intercept) 612.090  24.7405
>           Days         35.072   5.9221  0.066
>  Residual             654.941  25.5918
> # of obs: 180, groups: Subject, 18
> <snip>
>
>       I think we want the "Residual Variance" of 654.941 (Std.Dev. of
> 25.5918).  To get this, let's try "VarCorr":
>
> (vC.fm1 <- VarCorr(fm1))
> $Subject
> 2 x 2 Matrix of class "dpoMatrix"
>             (Intercept)     Days
> (Intercept)   612.09032  9.60428
> Days            9.60428 35.07165
>
> attr(,"sc")
> [1] 25.59182
>
>       Our desired 25.5918 is 'attr(,"sc")' here.  We get that as follows:
>
> (s. <- attr(vC.fm1, "sc"))
> [1] 25.59182
>
>       Next, by studying 'str(fm1)' and earlier emails you cite, we get
> the desired conditional posterior covariance matrix as follows:
>
> > (condPostCov <- (s.^2)*fm1@bVar$Subject)
> , , 1
>          [,1]       [,2]
> [1,] 145.7051 -21.444432
> [2,]   0.0000   5.312275
> <snip>
> , , 18
>          [,1]       [,2]
> [1,] 145.7051 -21.444432
> [2,]   0.0000   5.312275
>
>       This actually gives us only the upper triangular portion of the
> desired covariance matrices.  The following will make them all symmetric:
>
> condPostCov[2,1,] <- condPostCov[1,2,]
>
>       As a check, I computed the covariance matrix of the estimated
> random effects. To get this, I reviewed str(ranef(fm1)), which led me
> to the following:
>
>  var(ranef.fm1@.Data[[1]])
>             (Intercept)     Days
> (Intercept)   466.38503 31.04874
> Days           31.04874 29.75939
>
>       These numbers are all much larger than "condPostCov".  However,
> I believe this must be a random bounce -- unless it's a deficiency in
> my understanding (in which case, I hope someone will provide a
> correction).
>
>       Ulrich:  Would you mind checking this either with a published
> example or a Monte Carlo and reporting the results to us?
>
>       Viel Glück,
>       spencer graves
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Loading...