20190516_Biome_Mass_Multiphylo.Rmd

---
title: "20190324_habitat_mass_multiphylo"
author: "Nations"
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)

```

```{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
```