如何从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)

- b222

---

https://rdrr.io/github/JWiley/postMCMCglmm/man/ranef.MCMCglmm.html？ - Ben Bolker

我尝试使用 ranef(mc_mod2)，但出现了以下错误：Error in UseMethod("ranef") : no applicable method for 'ranef' applied to an object of class "MCMCglmm"。 - b222

你安装了那个包并加载它了吗？ - Ben Bolker

抱歉，我在那个网站上忽略了包名并假设它是mcmcglmm。非常感谢，Ben。 - b222

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

- Ben Bolker

---

这太棒了...谢谢，Ben。我的唯一问题是当我将随机效应指定为random = ~ idh(rem2):id + us(rem2):provider时，ranef函数似乎无法工作...我的意图是基于客户的种族-民族少数地位（白人/非白人）提取客户和提供者状态的随机效应。当我仅调用id和provider随机效应时，ranef函数运行良好，但是当我根据REM状态对它们进行分层时，ranef函数输出一个空矩阵。有什么想法吗？ - b222

最终目标是获得一种随机效果输出，看起来类似于：provider_1_white、provider_1_rem、provider_2_white、provider_2_rem 等等。也许这是不可能的吗？ - b222

我找到了一个解决方案......当随机效应分层（如上述模型）时，postMCMCglmm中的ranef()函数将无法使用（至少在此撰写时）。但是，可以通过使用str(mc_mod2)手动提取它们，然后使用cbind(B = posterior.mode(mc_mod2$Sol[,1:50]), CI = HPDinterval(mc_mod2$Sol[,1:50]))以可读格式呈现。 - b222

Answer 2

我忽略了需要安装的一个额外软件包（感谢Ben指出这一点）。

为了能够运行ranef()，只需安装postMCMCglmm软件包 - https://github.com/jwiley/postMCMCglmm/

#install.packages("devtools")
require(devtools)

install_github("JWiley/postMCMCglmm")

- b222