我想使用 mcmcsamp() 函数获得 lmer() 模型的方差分量的后验模拟。怎么做 ？

例如，下面是 lmer() 拟合的结果：

> 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

现在我运行 mcmcsamp() ：

> mm <- mcmcsamp(fit, n=15000) 

我希望剩余方差的模拟存储在“sigma”节点中，但这似乎不符合 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 

同样，我希望其他方差分量的模拟存储在“ST”节点中，但我得到了类似的观察结果。