我需要从模型中提取随机效应的后验估计和区间。
出于说明目的，与我正在使用的数据集类似的数据集将是 Base R 中的ChickWeight数据集。
我提取固定效应的后验估计和间隔的方式是这样的：
#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])
若要提取随机效应的估计值和 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。