r - 如何从 MCMCglmm 中提取随机效应？

Question

我正在寻找类似于 nlme、lme4 和 brms 中使用的 ranef() 的命令，它允许我在 MCMCglmm 模型中提取各个随机效应。在我的数据集中，我有 40 个提供商，我想提取每个提供商的随机效应并将它们绘制在毛毛虫图中。任何建议都会很棒。谢谢。

如果有帮助，这是我的 MCMCglmm 模型:

prior.3 <- list(R = list(R1 = list(V = diag(2), nu = 0.002)), 
                G = list(G1 = list(V = diag(2), nu = 0.002), 
                         G2 = list(V = diag(2), nu = 0.002))) 

mc_mod2 <- MCMCglmm(outcome ~ 1, data = filter(data, rem2 == "white" | rem2 == "rem"),
                  random = ~ idh(rem2):id + us(rem2):provider,
                  rcov = ~idh(rem2):units,
                  verbose = TRUE,
                  prior = prior.3,
                  family = "gaussian",
                  nitt = 100000, burnin = 5000,
                  pr = TRUE)

Answer 1

更详细一点，因为该包似乎没有内置毛毛虫图:请注意，在调用 MCMCglmm 时，您需要使用 pr=TRUE 以便存储随机效应值。

library(MCMCglmm)
data(PlodiaPO)  
model1 <- MCMCglmm(PO~1, random=~FSfamily, data=PlodiaPO, verbose=FALSE,
                 nitt=1300, burnin=300, thin=1,
                 pr=TRUE)
if (!require("postMCMCglmm")) {
    devtools::install_github("JWiley/postMCMCglmm")
    library("postMCMCglmm")
}

ranef() 似乎返回随机效应矩阵(行=级别，列=样本)。转换为具有均值和分位数的数据框:

qfun <- function(x,lev) unname(quantile(x,lev))
rsum <- as.data.frame(t(apply(ranef(model1),1,
      function(x) c(est=mean(x),
                    min=qfun(x,0.025),max=qfun(x,0.975)))))

绘图顺序:

rsum$term <- reorder(factor(rownames(rsum)),
                     rsum$est)

剧情:

library(ggplot2)
ggplot(rsum,aes(term,est))+
    geom_pointrange(aes(ymin=min,ymax=max))+
    coord_flip()