Wie extrahiere ich zufällige Effekte von MCMCglmm?

Question

Ich bin auf der Suche nach einem Befehl ähnlich wie ranef() in nlme, lme4 und brms verwendet, die mir erlauben, die einzelnen zufälligen Effekte in meinem MCMCglmm-Modell zu extrahieren. In meinem Datensatz habe ich 40 Anbieter und möchte die zufälligen Effekte für jeden Anbieter extrahieren und in einem Raupe-Plot darstellen. Irgendwelche Vorschläge wären großartig. Vielen Dank.

Falls es hilfreich ist, hier ist mein MCMCglmm Modell:

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

ich übersehen ein zusätzliches Paket, das installiert werden musste (Danke für diesen Hinweis, Ben).

zu können ranef() laufen, einfach das postMCMCglmm Paket installieren - https://github.com/jwiley/postMCMCglmm/

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

install_github("JWiley/postMCMCglmm") 

Answer 2

Etwas ausführlicher, da das Paket in haben gebaut Raupe Plots scheint nicht: beachten Sie pr=TRUE verwenden müssen beim Aufruf MCMCglmm, um die zufälligen Effekte Werte zu speichern.

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() erscheint eine Matrix der Zufallseffekten (rows = Ebenen, Säulen = Proben) zurückzukehren. Konvertieren in einen Datenrahmen mit mittlerem und Quantile:

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

Auftrag zum Plotten:

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

Grundstück:

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