I would like to get the posterior simulations of the variance components of a lmer() model with the mcmcsamp() function. How to do ?
For instance, below is the result of a lmer() fitting :
> fit Linear mixed model fit by REML Formula: y ~ 1 + (1 | Part) + (1 | Operator) + (1 | Part:Operator) Data: dat AIC BIC logLik deviance REMLdev 97.55 103.6 -43.78 89.18 87.55 Random effects: Groups Name Variance Std.Dev. Part:Operator (Intercept) 2.25724 1.50241 Part (Intercept) 3.30398 1.81769 Operator (Intercept) 0.00000 0.00000 Residual 0.42305 0.65043 Number of obs: 25, groups: Part:Operator, 15; Part, 5; Operator, 3
Now I run mcmcsamp() :
> mm <- mcmcsamp(fit, n=15000)
I expected that the simulations of the residual variance are stored in the "sigma" node but this does not seem to fit the results of lmer() :
> sigmasims <- mm@sigma[1,-(1:5000)] # discard first 5000 simulations (burn-in) > summary(sigmasims) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.8647 1.4960 1.7040 1.7460 1.9480 3.7920
Similarly I expected that the simulations of the other variance components are stored in the "ST" node but I get a similar observation.