提取 R 中 lme4 模型随机效应的后验估计和可信区间答案

【问题标题】：Extract posterior estimate and credible intervals for random effect for lme4 model in R提取 R 中 lme4 模型随机效应的后验估计和可信区间
【发布时间】：2020-03-20 04:17:42
【问题描述】：

我需要从我的模型中提取随机效应的后验估计和区间。

出于说明目的，与我正在使用的数据集类似的数据集是基础 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])

【问题讨论】：

str(sm.surv) 建议 sm.surv@ranef 和 sm.surv@fixef 是相当不同的
这是有道理的。在 r 中，fixef 是一个提取固定效应估计的函数，ranef 是一个提取随机效应估计的函数。
您可能希望使用merTools 包。或者，也许是 rstan 实现
为了扩展亨利的评论，str(smfixef.surv) 表明 smfixef.surv 是一个 100x5 矩阵。 str(ramranef.surv) 显示 smranef.surv 是一个包含 100x50x1 数组的列表。当我提取 smranef.surv=as.mcmc(smranef.surv[[1]][,,1]) 时，我能够申请 as.mcmc()。我认为这是可行的，因为 ?as.mcmc 显示该函数采用向量或矩阵，而不是列表或更高维数组。
@fausto.siegmund 这确实有助于它运行！但是，输出（可能是因为随机效应不同于固定效应）包括随机效应中每个水平的估计值和区间，而不是一个总体估计值和区间。关于如何正确合并这些值以便有一个估计值和一组随机效应间隔的任何想法？

标签： r lme4 confidence-interval mcmc credible-interval

【解决方案1】：

要提取随机效应的估计值和 95% CI，请使用以下代码：

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。

【讨论】：