0

我尝试使用 brms 包从系统发育模型中提取组级效应。

当我看到后验分布的平均值时，我设置了 95% 的区间以及哪些变量可以有效讨论。

model1tranformed <- ggs(bmodel)

exvars <- rownames(summary(bmodel)$fixed)
plots <- list()

for (i in 1:length(exvars)){
  p<-ggplot(filter(model1tranformed, Parameter == paste("b_", exvars[i], sep = "")),
            aes(x = value))+
    geom_density(fill  = "darkgreen", 
                 alpha = .5)+
    geom_vline(xintercept = 0, 
               col  = "darkred",
               size = 1)+
    scale_x_continuous(name   = "Value") + 
    geom_vline(xintercept = unlist(summary(bmodel)$fixed[i,3:4]),
               col = "blue",
               linetype = 2) +
    theme_light() +
    theme(text = element_text(size = 16))+
    labs(title = paste("pd_",exvars[i]))
  plots[[i]] = p
}

ggpubr::ggarrange(plotlist=plots,
                  ncol = 5, nrow = 4)

然后我想提取组级效果。我使用 ranef() 并可视化每个组的估计值。

qfun <- function(x,lev) unname(quantile(x,lev))
rsum<-list()
rsumm<-list()
plotr<-list()

for (i in 1:3){
  rsum[[i]] <- as.data.frame(t(apply(ranef(bmodel)[[1]][, , i],1,  
                                     function(x) c(est=mean(x),
                                                   min=qfun(x,0.025),max=qfun(x,0.975)))))
  
  rsum[[i]]$term <- reorder(factor(rownames(rsum[[i]])),
                            rsum[[i]]$est)
}

for (i in 1:3){
  rsumm[[i]]<-merge(rsum[[i]], dfpop_envsnu[,c(3,4,24)], by.x = "term", by.y = "sn", all.x = T)
  
  rsumm[[i]]<- rsumm[[i]] %>%
    distinct()
  rsumm[[i]]$japanese.name <- reorder(factor(rsumm[[i]]$japanese.name),
                           rsumm[[i]]$est)
}

for(i in 1:3){
  p <-ggplot(rsumm[[i]],aes(japanese.name,est))+
    geom_pointrange(aes(ymin=min,ymax=max))+
    coord_flip()+
    geom_hline(yintercept = 0,
               col  = "darkred",
               size = 1)+
    labs(title = paste(colnames(ranef(bmodel)[[1]][1 , , ])[i])) + 
    theme(plot.title = element_text(size=24),
          axis.title.y=element_blank(),
          axis.title.x=element_text(size=20))+
    ylab("Estimates")
  plotr[[i]] = p
}

ggpubr::ggarrange(plotlist=plotr,
                  ncol = 3, nrow = 1)

但我不知道哪个物种和变量有效或无效。是否有类似的方法来确认这些估计的 95% 区间？

4

0 回答 0