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