Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How do I extract random effects from MCMCglmm?

Tags:

r

bayesian

mcmc

I am looking for a command similar to ranef() used in nlme, lme4, and brms that will allow me to extract the individual random effects in my MCMCglmm model. In my dataset, I have 40 providers and I would like to extract the random effects for each provider and plot them in a caterpillar plot. Any suggestions would be great. Thank you.

In case it is helpful, here is my MCMCglmm model:

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)
like image 882
b222 Avatar asked Dec 01 '17 17:12

b222


1 Answers

A little more detail, since the package doesn't seem to have caterpillar plots built in: note you need to use pr=TRUE when calling MCMCglmm in order to store the random effects values.

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() appears to return a matrix of the random effects (rows=levels, columns=samples). Convert to a data frame with mean and quantiles:

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

Order for plotting:

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

Plot:

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

enter image description here

like image 168
Ben Bolker Avatar answered Sep 28 '22 16:09

Ben Bolker