Simulasi varian posterior dengan fungsi mcmcsamp

Saya ingin mendapatkan simulasi posterior komponen varians model lmer () dengan fungsi mcmcsamp (). Bagaimana melakukan ?

Misalnya, di bawah ini adalah hasil dari pas 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

Sekarang saya menjalankan mcmcsamp ():

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

Saya berharap bahwa simulasi varians residual disimpan dalam simpul "sigma" tetapi ini tampaknya tidak sesuai dengan hasil dari 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 

Demikian pula saya berharap bahwa simulasi komponen varians lainnya disimpan dalam "ST" node tetapi saya mendapatkan pengamatan serupa.

Jawaban singkat (ish) adalah itu

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

harus mengekstraksi rantai untuk efek tetap dan untuk efek acak dan varian residu dalam bentuk kerangka data.

Sebagai contoh:

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

Sayangnya vektor ini tidak mendapatkan nama yang berguna untuk komponen varians. Kamu bisa menggunakan

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

untuk memperbaiki hal ini (alih-alih mengkodekan posisi pada langkah terakhir Anda bisa melakukan sesuatu dengan -1:(length(fixef(fm2))))

Bagian lain dari jawaban ini adalah bahwa saya memiliki beberapa keraguan / pertanyaan serius tentang perilaku mcmcsamprantai, tetapi saya akan mencantumkan daftar: diskusi parsial / pendahuluan (dan mungkin salah!) Dari kebingungan saya diposting di http: //www.math.mcmaster.ca/~bolker/R/misc/mcmcsampex.pdf .