从MCMCglmm中提取随机效应

如何解决从MCMCglmm中提取随机效应

我有一个可以多次运行MCMCglmm的函数。

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

本节仅提取固定效果。

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

但是请注意，我不是在没有任何固定效果的情况下运行模型，因此输出为空。如何使用相同或相似的代码提取随机效果？

我想我可以将“解决方案”更改为其他内容，以从我运行的模型中提取随机效应而没有任何固定效应。

请注意，这是对上一个问题（带有示例df）的扩展-lapply instead of for loop for randomised hypothesis testing r

解决方法

使用broom.mixed::tidy是一个相对简单的方法。尚不清楚您是要提取顶级随机效果参数（即随机效果的方差）还是组级别效果的估算值的摘要。

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.

如果您想要组级效果，则需要effects="ran_vals"，但是必须使用pr=TRUE重新运行模型（或首先这样做），才能获得这些效果保存在模型对象中的效果：

 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