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

我正在寻找类似于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)

最佳答案

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

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()