Cleaning MCMCglmm model outputs

Dealing with multiple MCMCglmm model outputs?

Here is a a (reasonably) quick way of pulling out and collating model outputs from MCMCglmm objects. Those couple of functions allow you to grab the variable names, posterior means, CIs, effective sample sizes and pMCMC values for both fixed and random effects. A data.frame is returned and can then be used for tabulation or plotting purposes. It’s particularly handy when you have multiple models to deal with - that’s exactly why I needed it. It’s a quick and dirty way, but it worked great for what I needed it do it.

Code

Sort out your data and models

First, load your models and packages (you will need dplyr):

load("model1.RData")
load("model2.RData")
library(dplyr)  # need dplyr for the bind_rows function

Functions

The main function clean.MCMC() does most of the heavy lifting:

    clean.MCMC <- function(x) {
    sols <- summary(x)$solutions  ## pull out relevant info from model summary
    Gcovs <- summary(x)$Gcovariances
    Rcovs <- summary(x)$Rcovariances
    fixed <- data.frame(row.names(sols), sols, row.names = NULL)  ## convert to dataframes with the row.names as the first col
    random <- data.frame(row.names(Gcovs), Gcovs, row.names = NULL)
    residual <- data.frame(row.names(Rcovs), Rcovs, row.names = NULL)
    names(fixed)[names(fixed) == "row.names.sols."] <- "variable"  ## change the columns names to variable, so they all match
    names(random)[names(random) == "row.names.Gcovs."] <- "variable"
    names(residual)[names(residual) == "row.names.Rcovs."] <- "variable"
    fixed$effect <- "fixed"  ## add ID column for type of effect (fixed, random, residual)
    random$effect <- "random"
    residual$effect <- "residual"
    modelTerms <- as.data.frame(bind_rows(fixed, random, residual))  # merge it all together
}

The secondary function comes in when you have multiple models and need to add model names in as identifiers:

getName.MCMC <- function(x) deparse(substitute(x))  # add the model name

Try it!

Try it with one model:

oneModel <- clean.MCMC(model1)  # get all the info from summary(modelName)
oneModel$modelName <- getName.MCMC(model1)  # add the model's name in a new column
oneModel  # check out the created dataframe

And with multiple models… Note that the getName.MCMC won’t work nicely with a list unfortunately. There is probably a really easy way to fix that, but at the time a little workaround did the trick for me.

Manually create lists of input models:

dataList <- list(model1, model2)

Manually create lists of input model names:

dataListNames <- list("model1", "model2")

Get the clean.MCMC outputs and add modelName columns to each element for ID purposes:

readyList <- mapply(cbind, lapply(dataList, clean.MCMC), "modelName" = dataListNames, SIMPLIFY = F)

Turn the list of data.frames into one big data.frame:

mcmcOutputs <- as.data.frame(do.call(rbind, readyList), stringsAsFactors = FALSE)

Once you prepped the output as above you can use stargazer to create a nice summary table:

library(stargazer)
stargazer(mcmcOutputs, type = "text", summary = FALSE)

You can also plot lots of things, e.g. pull out the estimates for various models and plot those.