问题描述
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