Извлеките апостериорную оценку и вероятные интервалы для случайного эффекта для модели lme4 в R

Question

Мне нужно извлечь апостериорные оценки и интервалы для случайного эффекта из моей модели.

В иллюстративных целях набор данных, аналогичный тому, который я использую, будет ChickWeight набор данных в базе R.

Я извлекаю апостериорные оценки и интервалы для моих фиксированных эффектов следующим образом:

#load package
library(lme4)

#model
m.surv<-lmer(weight ~ Time + Diet + (1|Chick), data=ChickWeight)

#load packages
library(MCMCglmm)
library(arm)

#set up for fixed effects
sm.surv<-sim(m.surv)
smfixef.surv=sm.surv@fixef
smfixef.surv=as.mcmc(smfixef.surv)

#which gives
> posterior.mode(smfixef.surv)
(Intercept)        Time       Diet2  ... 
  8.5963329   8.7034260   5.1220436  ...
> HPDinterval(smfixef.surv)
                   lower      upper
(Intercept) -0.90309142 21.3617805
Time         8.42279728  9.0306337
Diet2       -6.84371527 35.1745980
...
attr(,"Probability")
[1] 0.95
>

Когда я пробую это для случайного эффекта (Chick), Я получаю следующую ошибку во второй строке кода:

smranef.surv=sm.surv@ranef
smranef.surv=as.mcmc(smranef.surv)

Ошибка в mcmc.list(x): аргументы должны быть объектами mcmc

Любые предложения о том, как я могу изменить свой код, чтобы извлечь эти значения для случайного эффекта?

Примечание для других пользователей: если модель была бы моделью MCMCglmm, значения апостериорной моды для выхода MCMC для случайных эффектов можно извлечь следующим образом:

posterior.mode(sm.surv$VCV[,1])
HPDinterval(sm.surv$VCV[,1])

Answer 1

Чтобы получить оценку и 95% доверительный интервал для случайных эффектов, вы используете следующий код:

sm.surv <-sim(m.surv)

#between Chick variance
bChick <-sm@ranef$Chick[,,1]
bvar<-as.vector(apply(bChick, 1, var)) #between ind variance posterior distribution
bvar<-as.mcmc(bvar)
posterior.mode(bvar) #mode of the distribution
HPDinterval(bvar)

Это даст вам:

>     posterior.mode(bvar)
     var1 
     501.24353 
>     HPDinterval(bvar)
      lower   upper
var1 412.36042 630.201
attr(,"Probability")
[1] 0.95

Это означает, что оценка составляет 501, нижний 95% интервал - 412, а верхний 95% интервал - 630.