# Posterior simulations of the variances with the mcmcsamp function

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.

-

as.data.frame(mm,type="varcov")


should extract the chains for the fixed effects and for the random-effect and residual variances in the form of a data frame.

For example:

library(lme4.0) ## I am using the r-forge version
fm2 <- lmer(Reaction ~ Days + (1|Subject) +
(0+Days|Subject), sleepstudy)
mm <- mcmcsamp(fm2,1000)
dd <- as.data.frame(mm,type="varcov")
burnin <- 100  ## probably unnecessary
summary(dd[-(1:burnin),])


Unfortunately this vector doesn't get useful names for the variance components. You can use

vnames <- c(names(getME(fm2,"theta")),"sigma^2")
names(dd)[3:5] <- vnames


to remedy this (instead of hard-coding the positions in the last step you could do something with -1:(length(fixef(fm2))))

The other part of this answer is that I am having some serious doubts/questions about the behavior of the mcmcsamp chains, but I will correspond off-list: a partial/preliminary (and possibly wrong!) discussion of my confusion is posted at http://www.math.mcmaster.ca/~bolker/R/misc/mcmcsampex.pdf .

-