extract random effects from MCMCglmm

Question

I have a function which runs an MCMCglmm a bunch of times.

shuffles <- 1:10
names(shuffles) <- paste0("shuffle_", shuffles)

library(MCMCglmm)
library(dplyr)
library(tibble)
library(purrr)

ddd <- purrr::map(shuffles,
                  ~ df %>%
                     mutate(Trait = sample(Trait)) %>%
                     MCMCglmm(fixed = Trait ~ 1,
                              random = ~ Year,
                              data = .,
                              family = "categorical",
                              verbose = FALSE)) %>%
   purrr::map( ~ tibble::as_tibble(summary(.x)$solutions, rownames = "model_term")) %>%
   dplyr::bind_rows(., .id = 'shuffle')
ddd

This section extracts fixed effects only.

(summary(.x)$Solutions, rownames = "model_term")

But note that I am not running a model without any fixed effects and so the output is empty. How can I extract random effects using the same or similar code?

I guess I can change 'solutions' to something else to extract random effects from a model I have run without any fixed effects.

Note that this is an extension to a previous question (with example df) here - lapply instead of for loop for randomised hypothesis testing r

This is an extension to a proviso question I asked (with example df) here - https://stackoverflow.com/questions/64287562/lapply-instead-of-for-loop-for-randomised-hypothesis-testing-r — Jamie Dunning, Oct 27 '20 at 20:45

score 2 · Answer 1 · answered Oct 27 '20 at 21:30

A relatively easy way to do this is with broom.mixed::tidy. It's not clear whether you mean you want to extract the summary for the top-level random effects parameters (i.e. the variances of the random effects), or for the estimates of the group-level effects.

library(broom.mixed)
tidy(m, effects="ran_pars")
##
##   effect   group    term                estimate std.error
## 1 ran_pars Year     var__(Intercept)     0.00212      629.
## 2 ran_pars Residual var__Observation 40465.         24211.

If you want the group-level effects you need effects="ran_vals", but you have to re-run your model with pr=TRUE (or do it that way in the first place) in order to have these effects saved in the model object:

 m <- MCMCglmm(Trait ~ ID, random = ~ Year, data = df, family = "categorical", pr=TRUE)
 tidy(m, effects="ran_vals")
 effect   group level term        estimate std.error
  <chr>    <chr> <chr> <chr>          <dbl>     <dbl>
1 ran_vals Year  1992  (Intercept)  2.65e-8      4.90
2 ran_vals Year  1993  (Intercept)  1.14e-8      6.23
3 ran_vals Year  1994  (Intercept)  1.28e-8      4.88
4 ran_vals Year  1995  (Intercept) -6.83e-9      5.31
5 ran_vals Year  1996  (Intercept) -1.36e-8      5.07
6 ran_vals Year  1997  (Intercept)  1.31e-8      5.24
7 ran_vals Year  1998  (Intercept) -2.80e-9      5.25
8 ran_vals Year  1999  (Intercept)  3.52e-8      5.68

Thanks Ben, but because Id' also like to run several permutations, I'd like to build this into the existing formula I have (linked on prev question). I am also getting an error Error: No tidy method for objects of class MCMCglmm — Jamie Dunning, Oct 28 '20 at 09:44

extract random effects from MCMCglmm

1 Answers1