CODEWITHSUNDEEP
rmcmc
Jamie Dunning
Jamie Dunning

Reputation: 153

extract random effects from MCMCglmm

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

Upvotes: 1

Views: 621

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226427

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

Upvotes: 2

Related Questions