# [R] bVar slot of lmer objects and standard errors

Ulrich Keller uhkeller at web.de
Thu Mar 2 11:25:48 CET 2006

```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 at 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.
>
> 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
>

```