Simulações posteriores das variações com a função mcmcsamp

Gostaria de obter as simulações posteriores dos componentes de variância de um modelo lmer () com a função mcmcsamp (). Como fazer ?

Por exemplo, abaixo está o resultado de um ajuste 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

Agora eu corro mcmcsamp ():

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

Eu esperava que as simulações da variação residual sejam armazenadas no nó "sigma", mas isso não parece se encaixar nos resultados de 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 

Da mesma forma, eu esperava que as simulações dos outros componentes de variância fossem armazenadas no nó "ST", mas recebo uma observação semelhante.

A resposta curta (ish) é que

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

deve extrair as cadeias para os efeitos fixos e para as variações residuais e de efeito aleatório na forma de um quadro de dados.

Por exemplo:

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),])

Infelizmente, esse vetor não obtém nomes úteis para os componentes de variação. Você pode usar

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

para remediar isso (em vez de codificar as posições na última etapa com a qual você poderia fazer algo -1:(length(fixef(fm2))))

A outra parte desta resposta é que estou tendo sérias dúvidas / dúvidas sobre o comportamento das mcmcsampcadeias, mas corresponderei fora da lista: uma discussão parcial / preliminar (e possivelmente errada!) Da minha confusão está publicada em http: //www.math.mcmaster.ca/~bolker/R/misc/mcmcsampex.pdf .