MCMCglmm: извлечение задней поверхности смоделированных фиксированных эффектов

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

В иллюстративных целях набор данных, аналогичный тому, который я использую, будет набором данных 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
>

Есть ли предложения о том, как я могу изменить свой код для извлечения SE для фиксированных эффектов (Time и Diet2)?


r bayesian mcmc standard-error
person Blundering Ecologist    schedule 16.12.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Вы можете извлечь их из сводной статистики smfixef.surv.

s <- summary(smfixef.surv)$statistics
s[grep("Time|Diet", rownames(s)), grep("SE", colnames(s))]
#         Naive SE Time-series SE
# Time  0.02046371     0.02046371
# Diet2 0.96785150     0.96785150
# Diet3 1.02553234     1.02553234
# Diet4 1.07393122     1.66512213
person jay.sf    schedule 16.12.2020
comment
В чем разница между наивным SE и SE временного ряда? - person Blundering Ecologist; 17.12.2020
comment
@BlunderingEcologist Согласно ?coda:::summary.mcmc, наивная стандартная ошибка среднего [игнорирует] автокорреляцию цепочки ... [тогда как] стандартная ошибка временного ряда [основывается] на оценке спектральной плотности на уровне 0. - person jay.sf; 17.12.2020