[R] 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 at 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 at .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 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