Chrotomyini_composite_072022.Rmd

---
title: "Chrotomyini_composite_072022"
author: "Stephanie M Smith"
date: '2022-07-12'
output: 
  html_document: 
    keep_md: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
pacman::p_load(install = F, "ape", "bayesplot","BiocManager", "brms", "broom", "devtools", "flextable","ggfortify", "ggrepel", "gtools", "lattice","loo", "nlme", "patchwork", "Rcpp", "remotes", "rstan", "smatr", "StanHeaders", "tidybayes", "viridis", "viridisLite", "pacman")
```

Data load in, phylogeny, transformations:
```{r}
d <- read.csv(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\05062022 Philippine Murids segmentation parameters and morphological data - TBA data total BoneJ (full).csv", header = T)

d <- d[d$tribe=="chroto",c(1:2, 4:23)]

d <- 
  d %>% 
  mutate(bvtv = as.numeric(bvtv))

d <- 
  d %>%
  mutate(mass_s = rethinking::standardize(log10(mass_g)),
         elev_s = rethinking::standardize(elev), 
         bvtv_s = rethinking::standardize(bvtv),
         tbth_s = rethinking::standardize(tbth),
         tbsp_s = rethinking::standardize(tbsp),
         conn_s = rethinking::standardize(conn),
         cond_s = rethinking::standardize(m_connd),
         cond_s2 = rethinking::standardize(connd), 
         da_s = rethinking::standardize(da))

# remove C. gonzalesi and R. isarogensis, singletons:
d <- 
  d %>% 
  filter(taxon!="Chrotomys_gonzalesi") %>% 
  filter(taxon!="Rhynchomys_isarogensis")

# Make categorical vars into factors
d <- 
  d %>%
  mutate(loco = factor(loco), 
         hab_simp = factor(hab_simp),
         genus = factor(genus))

# Means version
dmean <- d %>%
  group_by(taxon) %>% 
  summarize(mass_g = mean(mass_g),
            bvtv = mean(bvtv),
            tbth = mean(tbth),
            tbsp = mean(tbsp),
            m_connd = mean(m_connd))

# All log10 version
dlog <- d %>% 
  mutate(masslog = log10(mass_g),
         bvtvlog = log10(bvtv),
         tbthlog = log10(tbth),
         tbsplog = log10(tbsp),
         condlog = log10(m_connd)
  )
  

# Specify colors for genus-level plots:
cols = c("#00bbcd","#e6b531", "#adde79", "#2a098b","#ea908d")
# Order: Apomys, Archboldomys, Chrotomys, Rhynchomys, Soricomys

#Colors for Species-level plots:
cols2 = c("#00bbcd","#00bbcd","#00bbcd","#e6b531", "#adde79","#adde79","#adde79", "#2a098b","#ea908d", "#ea908d", "#ea908d")
# Order: A. banahao, A. datae, A. sierrae, A. maximus, C, mindorensis, C. silaceus, C. whiteheadi, S. kalinga, S. leonardocoi, S. montanus, R. labo

#### Phylogeny
ch.tre <- read.nexus(file = "G:\\My Drive\\Philippine rodents\\Chrotomys\\analysis\\SMS_PRUNED_and_COLLAPSED_03292022_OTUsrenamed_Rowsey_PhBgMCC_LzChrotomyini.nex")

ch <- ape::vcv.phylo(ch.tre, corr = T)

d <- 
  d %>% 
  mutate(phylo = taxon)

dlog <- 
  dlog %>% 
  mutate(phylo = taxon)
```

I. Bayesian phylogenetic multilevel models.
```{r}
# BV.TV, species
chr.1 <-
  brm(data = d, 
      family = student,
      bvtv_s ~ 0 + taxon + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.1")
print(chr.1)

# BV.TV, genus
chr.1.1 <-
  brm(data = d, 
      family = student,
      bvtv_s ~ 0 + genus + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.1.1")
print(chr.1.1)

chr.1.plot <- chr.1 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols2)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "BV.TV", y = "species")+
  ggtitle(label = "BV.TV by species")

chr.1.1.plot <- chr.1.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "BV.TV", y = "genus")+
  ggtitle(label = "BV.TV by genus")

chr.1.plot|chr.1.1.plot
```

```{r}
# Tb.Th, species
chr.2 <-
  brm(data = d, 
      family = student,
      tbth_s ~ 0 + taxon + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.2")
print(chr.2)

# Tb.Th, genus
chr.2.1 <-
  brm(data = d, 
      family = student,
      tbth_s ~ 0 + genus + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.2.1")
print(chr.2.1)

chr.2.plot <- chr.2 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols2)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "Tb.Th", y = "species")+
  ggtitle(label = "Tb.Th by species")

chr.2.1.plot <- chr.2.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "Tb.Th", y = "genus")+
  ggtitle(label = "Tb.Th by genus")

chr.2.plot|chr.2.1.plot

```

```{r}
# Tb.Sp, species
chr.3 <-
  brm(data = d, 
      family = student,
      tbsp_s ~ 0 + taxon + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.3")
print(chr.3)

# Tb.Sp, genus
chr.3.1 <-
  brm(data = d, 
      family = student,
      tbsp_s ~ 0 + genus + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.3.1")
print(chr.3.1)

chr.3.plot <- chr.3 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols2)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "Tb.Sp", y = "species")+
  ggtitle(label = "Tb.Sp by species")

chr.3.1.plot <- chr.3.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "Tb.Sp", y = "genus")+
  ggtitle(label = "Tb.Sp by genus")

chr.3.plot|chr.3.1.plot
```

```{r}
# conn.D, species
chr.4 <-
  brm(data = d, 
      family = student,
      cond_s ~ 0 + taxon + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.4")
print(chr.4)


# conn.D, genus
chr.4.1 <-
  brm(data = d, 
      family = student,
      cond_s ~ 0 + genus + mass_s + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.95), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.4.1")
print(chr.4.1)


chr.4.plot <- chr.4 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .89) +
  scale_fill_manual(values = cols2)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "Conn.D", y = "species")+
  ggtitle(label = "Conn.D by species")

chr.4.1.plot <- chr.4.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye(aes(fill = .variable), 
               point_fill = "#000000", 
               shape = 21, 
               point_size = 3, 
               point_color = "#FFFFFF",
               interval_size = 10,
               interval_color = "grey40",
               .width = .95) +
  scale_fill_manual(values = cols)+
  geom_vline(xintercept = 0, linetype = "dashed")+
  theme(legend.position = "none", 
        plot.title = element_text(size = 9))+
  labs(x = "Conn.D", y = "genus")+
  ggtitle(label = "Conn.D by genus")

chr.4.plot|chr.4.1.plot
```

II. Comparison of TB metrics between taxa (substrate use proxy).

```{r}
#### BV.TV ####
bvtv.comp.s <-chr.1 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-species BV.TV difference distributions")+
  labs(x = "difference", y = "comparison")

bvtv.comp.g <-chr.1.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-genus BV.TV difference distributions")+
  labs(x = "difference", y = "comparison")
bvtv.comp.g
bvtv.comp.s
```

```{r}
#### Tb.Th ####
tbth.comp.s <-chr.2 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-species Tb.Th difference distributions")+
  labs(x = "difference", y = "comparison")

tbth.comp.g <-chr.2.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-genus Tb.Th difference distributions")+
  labs(x = "difference", y = "comparison")
tbth.comp.g
tbth.comp.s
```

```{r}
#### Tb.Sp ####
tbsp.comp.s <-chr.3 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-species Tb.Sp difference distributions")+
  labs(x = "difference", y = "comparison")

tbsp.comp.g <-chr.3.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-genus Tb.Sp difference distributions")+
  labs(x = "difference", y = "comparison")
tbsp.comp.g
tbsp.comp.s
```

```{r}
#### Conn.D ####
cond.comp.s <-chr.4 %>%
  gather_draws(b_taxonApomys_banahao, b_taxonApomys_datae, b_taxonApomys_sierrae, b_taxonArchboldomys_maximus, b_taxonChrotomys_mindorensis, b_taxonChrotomys_silaceus, b_taxonChrotomys_whiteheadi, b_taxonRhynchomys_labo, b_taxonSoricomys_kalinga, b_taxonSoricomys_leonardocoi, b_taxonSoricomys_montanus) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_taxon", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-species Conn.D difference distributions")+
  labs(x = "difference", y = "comparison")

cond.comp.g <-chr.4.1 %>%
  gather_draws(b_genusApomys,b_genusArchboldomys, b_genusChrotomys, b_genusRhynchomys, b_genusSoricomys) %>%
  compare_levels(.value, by = .variable) %>%
  mutate(.variable = str_replace_all(.variable, "b_genus", "")) %>%
  ungroup() %>%
  mutate(loco = reorder(.variable, .value)) %>%
  ggplot(aes(y = .variable, x = .value)) +
  stat_halfeye() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ggtitle(label = "Between-genus Conn.D difference distributions")+
  labs(x = "difference", y = "comparison")
cond.comp.g
cond.comp.s
```

III. Allometry with brms, smatr, and pgls.

```{r}
#### BV.TV ####
# brms
chr.5 <-
  brm(data = dlog, 
      family = student,
      bvtvlog ~ 1 + masslog + (1|gr(phylo, cov = ch)),
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.5")
print(chr.5)

# smatr
bvtv.sma <- sma(bvtv~mass_g, data = dmean, slope.test = 0, log = "xy", robust = T)
summary(bvtv.sma)

# PGLS
bvtv.pgls <-  gls(log10(bvtv) ~ log10(mass_g), correlation = corBrownian(phy = ch.tre, form = ~taxon), data = d, method = "ML")
coef(bvtv.pgls)

#### Tb.Th ####
# brms
chr.6 <-
  brm(data = dlog, 
      family = student,
      tbthlog ~ 1 + masslog + (1|gr(phylo, cov = ch)),
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.6")
print(chr.6)

# smatr
tbth.sma <- sma(tbth~mass_g, data = dmean, slope.test = 0, log = "xy", robust = T)
summary(tbth.sma)

# PGLS
tbth.pgls <-  gls(log10(tbth) ~ log10(mass_g), correlation = corBrownian(phy = ch.tre, form = ~taxon), data = d, method = "ML")
coef(tbth.pgls)

#### Tb.Sp ####
# brms
chr.7 <-
  brm(data = dlog, 
      family = student,
      tbsplog ~ 1 + masslog + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.98), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.7")
print(chr.7)

# smatr
tbsp.sma <- sma(tbsp~mass_g, data = dmean, slope.test = 0, log = "xy", robust = T)
summary(tbsp.sma)

# PGLS
tbsp.pgls <-  gls(log10(tbsp) ~ log10(mass_g), correlation = corBrownian(phy = ch.tre, form = ~taxon), data = d, method = "ML")
coef(tbsp.pgls)

#### conn.D ####
# brms
chr.8 <-
  brm(data = dlog, 
      family = student,
      condlog ~ 1 + masslog + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.9), #inserted to decrease the number of divergent transitions here
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = b),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.8")
print(chr.8)

# smatr
cond.sma <- sma(m_connd~mass_g, data = dmean, slope.test = 0, log = "xy", robust = T)
summary(cond.sma)

# PGLS
cond.pgls <-  gls(log10(m_connd) ~ log10(mass_g), correlation = corBrownian(phy = ch.tre, form = ~taxon), data = d, method = "ML")
coef(cond.pgls)
```


IV. Phylogenetic signal of trabecular bone metrics.

```{r}
# Establish hypothesis to be tested
hyp <- "sd_phylo__Intercept^2 / (sd_phylo__Intercept^2 + sigma^2) = 0"

# BV.TV
ps.bvtv.1 <- hypothesis(chr.1.1, hyp, class = NULL)
ps.bvtv.1 # genus
ps.bvtv <- hypothesis(chr.1, hyp, class = NULL)
ps.bvtv # species

# Tb.Th
ps.tbth.1 <- hypothesis(chr.2.1, hyp, class = NULL)
ps.tbth.1 # genus
ps.tbth <- hypothesis(chr.2, hyp, class = NULL)
ps.tbth # species

# Tb.Sp
ps.tbsp.1 <- hypothesis(chr.3.1, hyp, class = NULL)
ps.tbsp.1 # genus
ps.tbsp <- hypothesis(chr.3, hyp, class = NULL)
ps.tbsp # species

# conn.D
ps.cond.1 <- hypothesis(chr.4.1, hyp, class = NULL)
ps.cond.1  # genus
ps.cond <- hypothesis(chr.4, hyp, class = NULL)
ps.cond # species

#### Plots: species-level ####

bvtv.stpl <- ggplot() +
  geom_density(aes(x = ps.bvtv$samples$H1), fill = "#fc2187", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: bone volume fraction")

tbth.stpl <- ggplot() +
  geom_density(aes(x = ps.tbth$samples$H1), fill = "#66116a", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: trabecular thickness")

tbsp.stpl <- ggplot() +
  geom_density(aes(x = ps.tbsp$samples$H1), fill = "#9f9ad2", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: trabecular separation")

cond.stpl <- ggplot() +
  geom_density(aes(x = ps.cond$samples$H1), fill = "#ecd156", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: connectivity density")
bvtv.stpl/cond.stpl|tbth.stpl/tbsp.stpl
```

```{r}
#### Plots: genus-level ####
bvtv.stpl.1 <- ggplot() +
  geom_density(aes(x = ps.bvtv.1$samples$H1), fill = "#fc2187", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: bone volume fraction")

tbth.stpl.1 <- ggplot() +
  geom_density(aes(x = ps.tbth.1$samples$H1), fill = "#66116a", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: trabecular thickness")

tbsp.stpl.1 <- ggplot() +
  geom_density(aes(x = ps.tbsp.1$samples$H1), fill = "#9f9ad2", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: trabecular separation")

cond.stpl.1 <- ggplot() +
  geom_density(aes(x = ps.cond.1$samples$H1), fill = "#ecd156", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: connectivity density")
bvtv.stpl.1/cond.stpl.1|tbth.stpl.1/tbsp.stpl.1
```


V. Connectivity density (Conn.D) vs BV.TV
```{r}
cor.test(d$cond_s,d$bvtv_s) # cor = -0.16, p = 0.18. CIs = -0.38-0.080. Weak.
ggplot(aes(y = cond_s, x = bvtv_s, color = genus), data = d)+
  scale_fill_manual(values = cols)+
  geom_point()
  
```

VI. Phylogenetic signal of mass alone.

```{r}
chr.9 <-   
  brm(data = d, 
      family = student,
      mass_s ~ 0 + (1|gr(phylo, cov = ch)),
      control = list(adapt_delta = 0.98), #inserted to decrease the number of divergent transitions
      prior = c(
        prior(gamma(2, 0.1), class = nu),
        prior(normal(0, 1), class = sd),
        prior(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\chr.9")
print(chr.9)
```

```{r}
ps.mass <- hypothesis(chr.9, hyp, class = NULL)
ps.mass #lambda = 0.95, CIs: 0.91-0.98
mass.stpl <- ggplot() +
  geom_density(aes(x =ps.mass$samples$H1), fill = "#fc2187", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: mass")

mass.stpl
```