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([EMAIL PROTECTED],,]*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)[EMAIL PROTECTED])
> , , 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([EMAIL PROTECTED])
> (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
>
______________________________________________
R-help@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
