--- title: "20190324_habitat_mass_multiphylo" author: "" date: "5/16/2019" output: html_document --- ### This script estimates the phylogenetically corrected body mass for 3 different Australia biomes. ####This uses 100 randomly sampled phylogenetic trees from a posterior distribution generate by Smissen and Rowe, 2018 #### Warning, takes a long time to run! Set working directory ```{r setup, include=FALSE} knitr::opts_knit$set(root.dir = "~/User....") ``` required packages ```{r, message=FALSE} library(tidyverse) library(brms) library(MCMCglmm) library(tidybayes) library(bayesplot) library(phytools) getwd() ``` Acquire and mean center data load tree and trees, and subsample ```{r} tree <- read.tree("Pseudomys_Division_Tree.nwk") trees <- read.nexus("5000_Posterior.trees") tree_samp <- sample(trees, size = 100) #tree <- read.tree("") ``` ```{r} dat <- read.csv("Pseudomys_Division_data.csv", header=TRUE, stringsAsFactors=FALSE) dat$log_mass <-log(dat$mass) dat$cmass <- dat$mass - mean(dat$mass) dat$temp <- dat$Temperature_annual_mean - mean(dat$Temperature_annual_mean) dat$ar <- dat$Aridity_index_annual_mean - mean(dat$Aridity_index_annual_mean) dat$pam <- dat$Precip_annual_mean - mean(dat$Precip_annual_mean) dat$el <- dat$Elevation - mean(dat$Elevation) dat$pas <- dat$Precip_annual_seasonality_ratio - mean(dat$Precip_annual_seasonality_ratio) dat$pwm <- dat$Precip_wettest_month - mean(dat$Precip_wettest_month) dat$tar <- dat$Temperature_annual_range - mean(dat$Temperature_annual_range) dat <- dat[order(factor(dat$Species, levels = tree$tip.label)), ] ``` For 1 tree ```{r} inv.phylo <- MCMCglmm::inverseA(tree, nodes = "TIPS", scale = TRUE) A <- solve(inv.phylo$Ainv) rownames(A) <- rownames(inv.phylo$Ainv) ``` Generate a matrix of 100 covariance matrices ```{r} inverse<-list() As <- list() for (i in 1:length(tree_samp) ) { inverse[[i]] <- inverseA(tree_samp[[i]], nodes = "TIPS", scale = TRUE) x <- solve(inverse[[i]]$Ainv) rownames(x) <- rownames(inverse[[i]]$Ainv) As[[i]] <- x } ``` ### Multiphylo test ```{r} h.fit <- brm( log_mass ~ 0 + hab_name + (1 | phylo), data = dat, family = student(link = "identity"), cov_ranef = list(phylo = As[[1]]), iter = 5000, save_all_pars = TRUE, prior = c( set_prior("normal(0,1.5)", class = "b"), set_prior("normal(0,1.5)", class = "sd")), control = list(adapt_delta = 0.99, max_treedepth = 15), cores = parallel::detectCores() ) h.fits <- vector("list", 100) for (i in seq_along(h.fits)) { h.fits[[i]] <- update(h.fit, cov_ranef = list(phylo = As[[i]], cores = parallel::detectCores()) ) } ``` #### Insert file name here to save Rds ```{r} h.fits_comb <- combine_models(h.fits[[i]], mlist = h.fits) saveRDS(h.fits_comb, file = "~/User...../file.Rds", compress = F) ## To read in the rds file: # brms_output <- readRDS(file = "brms_saved_file.Rds") summary(h.fits_comb) plot(h.fits_comb, ask = F) ``` ### Plotting sposterior distributions ```{r} A.Me <- posterior_samples(h.fits_comb) %>% transmute(dif1 = b_hab_nameArid - b_hab_nameMesic) %>% ggplot(aes(x = dif1, y = 0)) + geom_halfeyeh(fill = "firebrick4", point_interval = median_qi, .width = .95) + scale_y_continuous(NULL, breaks = NULL) + labs(subtitle = "Difference score, Arid - Mesic", x = expression("Arid - Mesic")) + theme_bw() + theme(panel.grid = element_blank()) A.Mo <- posterior_samples(h.fits_comb) %>% transmute(dif2 = b_hab_nameArid - b_hab_nameMonsoon) %>% ggplot(aes(x = dif2, y = 0)) + geom_halfeyeh(fill = "firebrick4", point_interval = median_qi, .width = .95) + scale_y_continuous(NULL, breaks = NULL) + labs(subtitle = "Difference score, Arid - Monsoon", x = expression("Arid - Monsoon")) + theme_bw() + theme(panel.grid = element_blank()) Mo.Me <- posterior_samples(h.fits_comb) %>% transmute(dif3 = b_hab_nameMonsoon - b_hab_nameMesic) %>% ggplot(aes(x = dif3, y = 0)) + geom_halfeyeh(fill = "firebrick4", point_interval = median_qi, .width = .95) + scale_y_continuous(NULL, breaks = NULL) + labs(subtitle = "Difference score, Monsoon - Mesic", x = expression("Monsoon - Mesic")) + theme_bw() + theme(panel.grid = element_blank()) ``` Plotting Pagel's Lambda value ```{r} hyp <- paste( "sd_phylo__Intercept^2 /", "(sd_phylo__Intercept^2 + sigma^2) = 0" ) (hyp <- hypothesis(h.fits_comb, hyp, class = NULL)) hyppp <- plot(hyp, plot = F, theme = theme_get())[[1]] hyp_1 <- hyp$samples which.max(density(hyp_1$H1)$y) density(hyp_1$H1)$x[500] ``` ```{r} hyp_1p <- ggplot(hyp_1, aes(H1, y=0)) + geom_density(color = "black", fill = "firebrick4") + geom_vline(xintercept = density(hyp_1$H1)$x[500], size = 0.3) + geom_vline(xintercept = mean(hyp_1$H1), linetype = "dashed", size = 0.3) + theme_bw() + labs(x = expression("Body Mass "*lambda), subtitle = "Difference score, Arid - Monsoon",center = TRUE) + ylim(0,40) + xlim(0,1.1) p_plot <- plot_grid(A.Me, A.Mo, Mo.Me, hyppp, nrow = 2, ncol = 2) p_plot ```