Simulaciones posteriores de las varianzas con la función mcmcsamp

Question

Me gustaría obtener las simulaciones posteriores de los componentes de varianza de un modelo lmer() con la función mcmcsamp(). ¿Cómo hacerlo?

Por ejemplo, a continuación se muestra el resultado de un ajuste de 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

Ahora ejecuto mcmcsamp() :

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

Esperaba que las simulaciones de la varianza residual se almacenaran en el nodo "sigma", pero esto no parece ajustarse a los 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 

Del mismo modo, esperaba que las simulaciones de los demás componentes de la varianza se almacenaran en el nodo "ST", pero obtengo una observación similar.

Answer 1

La respuesta corta es que

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

debe extraer las cadenas para los efectos fijos y para las varianzas de efectos aleatorios y residuales en forma de marco de datos.

Por ejemplo:

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

Por desgracia, este vector no tiene nombres útiles para los componentes de la varianza. Puede utilizar

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

para remediarlo (en lugar de codificar las posiciones en el último paso, podría hacer algo con -1:(length(fixef(fm2))) )

La otra parte de esta respuesta es que estoy teniendo algunas serias dudas / preguntas sobre el comportamiento de la mcmcsamp cadenas, pero corresponderé fuera de la lista: una discusión parcial/preliminar (¡y posiblemente errónea!) de mi confusión está publicada en http://www.math.mcmaster.ca/~bolker/R/misc/mcmcsampex.pdf .