Chrotomyini_sandbox_analyses_06032022.Rmd

---
title: 'Chrotomys analysis: Student priors, genus/sp/etc'
author: "S.M. Smith"
date: "5/2/2022"
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", "dagitty", "devtools", "flextable", "ggdark", "ggfortify", "ggmcmc", "ggrepel", "gtools", "lattice","loo", "patchwork", "rcartocolor", "Rcpp", "remotes", "rstan", "StanHeaders", "statebins", "tidybayes", "viridis", "viridisLite", "pacman")
```

Load up Chrotomyini trabecular bone architecture (TBA) data and standardize variables:
```{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))

# Specify colors for plots:
cols = c("#86acca","#ab685b", "#3370a3", "#1f314c","#5a9fa8")
```

Load in phylogeny: 
REMEMBER: A <- ape::vcv.phylo(phylo), add corr = T if your tree is NOT SCALED TO 1. 
```{r}
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)
```

From existing fits: call up some models for comparison. These are all BV.TV.
```{r}

# Just mass, wide gamma for nu
ch.70 <- 
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.70")
ch.70 <- add_criterion(ch.70, c("loo", "waic"))


# Just genus, wide gamma for nu
ch.71 <- 
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.71")
ch.71 <- add_criterion(ch.71, c("loo", "waic"))

# Just genus, nu fixed at 2
ch.71b <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.71b")
ch.71b <- add_criterion(ch.71b, c("loo", "waic"))

# By genus and mass
ch.72 <- 
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.72")
ch.72 <- add_criterion(ch.72, c("loo", "waic"))

# By genus and mass with fixed nu
ch.72b <- 
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.72b")
ch.72b <- add_criterion(ch.72b, c("loo", "waic"))

# By genus, mass, and BOTH phylogeny/intraspecific variance
ch.74 <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.74")
ch.74 <- add_criterion(ch.74, c("loo", "waic"))

# By genus, mass, and phylogeny only
ch.74.1 <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.74.1")
mcmc_plot(ch.74.1, pars = "^b_")
ch.74.1 <- add_criterion(ch.74.1, c("loo", "waic"))

# By genus, mass, and intraspecific variance only
ch.74.2 <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.74.2")
ch.74.2 <- add_criterion(ch.74.2, c("loo", "waic"))

# by SPECIES, mass, and phylo
ch.74.3 <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.74.3")
print(ch.74.3)
mcmc_plot(ch.74.3, pars = "^b_")
ch.74.3 <- add_criterion(ch.74.3, c("loo", "waic"))


# By genus, mass, phylogeny/intraspcific var, fixed nu = 2
ch.74b <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.74b")
ch.74b <- add_criterion(ch.74b, c("loo", "waic"))

# By genus, mass, phylogeny/intraspcific var, with estimated nu, more robust gamma(4,1) as demonstrated by Kurz - a slightly more restrictive prior but not as much as the fixed one:

ch.74c <-
  brm(file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.74c")
ch.74c <- add_criterion(ch.74c, c("loo", "waic"))

b_estimates <-
  tibble(model = c("ch.74", "ch.74b", "ch.74c")) %>% 
  mutate(fit = map(model, get)) %>% 
  mutate(posterior_summary = map(fit, ~posterior_summary(.) %>% 
                                   data.frame() %>% 
                                   rownames_to_column("term"))) %>% 
  unnest(posterior_summary) %>% 
  select(-fit) %>% 
  arrange(term)

```

Ok try some other things. Model just based on phylogeny:
```{r}
ch.75 <- 
  brm(data = d, 
      family = student,
      bf(bvtv_s ~ 1+ (1|gr(phylo, cov = ch)), nu = 2),
      prior = c(
        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\\ch.75")

ch.75 <- add_criterion(ch.75, c("loo", "waic"))

loo_compare(ch.75, ch.74, ch.74.1, ch.74.2, ch.74.3, ch.70, ch.71, ch.71b)%>% 
  print(simplify = F)
```

Try some more species things: 
```{r}

# by species tbth
ch.76 <-
  brm(data = d, 
      family = student,
      tbth_s ~ 0 + taxon + mass_s + (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\\ch.76")
ch.76 <- add_criterion(ch.76, c("loo", "waic"))

# by phylo only tbth
ch.76.1 <- 
  brm(data = d, 
      family = student,
      bf(tbth_s ~ 1+ (1|gr(phylo, cov = ch)), nu = 2),
      prior = c(
        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\\ch.76.1")
ch.76.1 <- add_criterion(ch.76.1, c("loo", "waic"))

hyp <- "sd_phylo__Intercept^2 / (sd_phylo__Intercept^2 + sigma^2) = 0"
h.bvtv <- hypothesis(ch.76, hyp, class = NULL)
h.bvtv

# By mass only
ch.76.2 <- 
  brm(data = d, 
      family = student,
      tbth_s ~ 1+mass_s,
      prior = c(prior(normal(0, 1), class = b),
                prior(gamma(2, 0.1), class = nu),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 5,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.76.2")

ch.76.2 <- add_criterion(ch.76.2, c("loo", "waic"))

loo_compare(ch.76, ch.76.1, ch.76.2) %>% 
  print(simplify = F)
print(ch.76)
mcmc_plot(ch.76, pars = "^b_")
```

```{r}
# by species tbsp
ch.77 <-
  brm(data = d, 
      family = student,
      tbsp_s ~ 0 + taxon + mass_s + (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\\ch.77")
print(ch.77)
mcmc_plot(ch.77, pars = "^b_")



# by phylo only tbsp
ch.77.1 <- 
  brm(data = d, 
      family = student,
      bf(tbsp_s ~ 1+ (1|gr(phylo, cov = ch)), nu = 2),
      prior = c(
        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\\ch.77.1")
ch.77.1 <- add_criterion(ch.77.1, c("loo", "waic"))

# By mass only
ch.77.2 <- 
  brm(data = d, 
      family = student,
      tbsp_s ~ 1+mass_s,
      prior = c(prior(normal(0, 1), class = b),
                prior(gamma(2, 0.1), class = nu),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 5,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.77.2")

ch.77.2 <- add_criterion(ch.77.2, c("loo", "waic"))

# By genus tbsp
ch.77.3 <-
  brm(data = d, 
      family = student,
      tbsp_s ~ 0 + genus + mass_s + (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\\ch.77.3")
print(ch.77.3)
hyp <- "sd_phylo__Intercept^2 / (sd_phylo__Intercept^2 + sigma^2) = 0"
h.bvtv <- hypothesis(ch.77.3, hyp, class = NULL)
h.bvtv
ch.77.3 <- add_criterion(ch.77.3, c("loo", "waic"))
plot(ch.77.3)
ph.bvtv.pl <- ggplot() +
  geom_density(aes(x = h.bvtv$samples$H1), fill = "red", alpha = 0.5) +
  theme_bw() +
  xlim(0,1) +
  labs(y = "density", x = "lambda: Trabecular spacing")
ph.bvtv.pl

loo_compare(ch.77, ch.77.1, ch.77.2, ch.77.3) %>% 
  print(simplify = F)

```

```{r}
# by species cond
ch.78 <-
  brm(data = d, 
      family = student,
      cond_s ~ 0 + taxon + mass_s + (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\\ch.78")
print(ch.78)
mcmc_plot(ch.78, pars = "^b_")
ch.78 <- add_criterion(ch.78, c("loo", "waic"))


# by phylo only cond
ch.78.1 <- 
  brm(data = d, 
      family = student,
      bf(cond_s ~ 1+ (1|gr(phylo, cov = ch)), nu = 2),
      prior = c(
        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\\ch.78.1")
ch.78.1 <- add_criterion(ch.78.1, c("loo", "waic"))
print(ch.78.1)

# By mass only
ch.78.2 <- 
  brm(data = d, 
      family = student,
      cond_s ~ 1+mass_s,
      prior = c(prior(normal(0, 1), class = b),
                prior(gamma(2, 0.1), class = nu),
                prior(exponential(1), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 5,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.78.2")

ch.78.2 <- add_criterion(ch.78.2, c("loo", "waic"))

# By genus tbsp
ch.78.3 <-
  brm(data = d, 
      family = student,
      cond_s ~ 0 + genus + mass_s + (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\\ch.78.3")
print(ch.78.3)

ch.78.3 <- add_criterion(ch.78.3, c("loo", "waic"))

loo_compare(ch.78, ch.78.1, ch.78.2, ch.78.3) %>% 
  print(simplify = F)

```

Some models by species and mass only to plot on the phylogeny as a continuous trait. Is this allowed? I don't know for sure but I'm gonna try it.
```{r}
# by species tbth no phy
ch.75.4 <-
  brm(data = d, 
      family = student,
      bvtv_s ~ 0 + taxon + mass_s,
      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(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.75.4")
print(ch.75.4)
mcmc_plot(ch.75.4, pars = "^b_")
ch.75.4 <- add_criterion(ch.75.4, c("loo", "waic"))

# by species tbth no phy
ch.76.4 <-
  brm(data = d, 
      family = student,
      tbth_s ~ 0 + taxon + mass_s,
      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(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.76.4")
print(ch.76.4)
mcmc_plot(ch.76.4, pars = "^b_")
ch.76.4 <- add_criterion(ch.76.4, c("loo", "waic"))

# by species tbsp no phy
ch.77.4 <-
  brm(data = d, 
      family = student,
      tbsp_s ~ 0 + taxon + mass_s,
      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(exponential(1), class = sigma)
        ),
      data2 = list(ch = ch),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.77.4")
print(ch.77.4)
mcmc_plot(ch.77.4, pars = "^b_")
ch.77.4 <- add_criterion(ch.77.4, c("loo", "waic"))

# by species cond no phy
ch.78.4 <-
  brm(data = d, 
      family = student,
      cond_s ~ 0 + taxon + mass_s,
      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(exponential(1), class = sigma)
        ),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      file = "G:\\My Drive\\Philippine rodents\\chrotomyini\\fits\\ch.78.4")
print(ch.78)
mcmc_plot(ch.78.4, pars = "^b_")
ch.78 <- add_criterion(ch.78, c("loo", "waic"))
```

How about some PCAs? Just real quick.

```{r}
d_pc <- d %>% 
  select(c("taxon", "specno","genus", (ends_with("_s")))) %>% 
  select(!c("elev_s", "da_s", "conn_s")) %>% 
  column_to_rownames("specno")
pc_s <- prcomp(d_pc[,3:7])

pc12 <- autoplot(pc_s, x = 1, y = 2, data = d_pc, colour = "genus", size = 3, loadings.label = T, loadings = T, loadings.colour = "#000000", loadings.label.colour = "#000000",scale = 0)

pc12species <- autoplot(pc_s, x = 1, y = 2, data = d_pc,  colour = "taxon", size = 3, loadings.label = T, loadings = T, shape = "genus", scale = 0)

pc32 <- autoplot(pc_s, x = 3, y = 2, data = d_pc, colour = "genus", size = 3, loadings.label = T, loadings = T, loadings.colour = "#000000", loadings.label.colour = "#000000", scale = 0)

pc12|pc32

pcall <- pc_s$x %>% 
  as.data.frame() %>% 
  mutate(genus = d_pc$genus,
         taxon = d_pc$taxon)

pccolsplot <- autoplot(mapping = aes(lwd = 10), pc_s, x = 1, y = 2, data = d_pc, shape = F, label =  F, loadings.label = T, loadings = T, loadings.colour = "#000000", loadings.label.colour = "#000000", scale = 0) +
  geom_point(aes(PC1, y=PC2, color = taxon),  data = pcall, size = 8)+
  scale_color_manual(values = cols2)
pccolsplot



```

Can I plot JUST the Chrotomys guys?

```{r}
pcchrot <- pc_s$x %>% 
  as.data.frame() %>% 
  mutate(genus = d_pc$genus,
         taxon = d_pc$taxon) %>% 
  filter(genus=="Chrotomys")

pcchplot <- autoplot(pc_s, x = 1, y = 2, data = d_pc, shape = F, label = F, loadings.label = T, loadings = T, loadings.colour = "#000000", loadings.label.colour = "#000000", scale = 0) +
  geom_point(aes(PC1, y=PC2, color = taxon), data = pcchrot, size = 3)

pcchplot
```

May as well do ones for Apomys and Soricomys. 
Apomys:
```{r}
pcap <- pc_s$x %>% 
  as.data.frame() %>% 
  mutate(genus = d_pc$genus,
         taxon = d_pc$taxon) %>% 
  filter(genus=="Apomys")

pcapplot <- autoplot(pc_s, x = 1, y = 2, data = d_pc, shape = F, label = F, loadings.label = T, loadings = T, loadings.colour = "#000000", loadings.label.colour = "#000000", scale = 0) +
  geom_point(aes(PC1, y=PC2, color = taxon), data = pcap, size = 3)

pcapplot
```

Soricomys:
```{r}
pcsor <- pc_s$x %>% 
  as.data.frame() %>% 
  mutate(genus = d_pc$genus,
         taxon = d_pc$taxon) %>% 
  filter(genus=="Soricomys")

pcsorplot <- autoplot(pc_s, x = 1, y = 2, data = d_pc, shape = F, label = F, loadings.label = T, loadings = T, loadings.colour = "#000000", loadings.label.colour = "#000000", scale = 0) +
  geom_point(aes(PC1, y=PC2, color = taxon), data = pcsor, size = 3)

pcsorplot
```

How about some allometry in case someone asks?
```{r}
ggplot(aes(x=bvtv, y = mass_g), data = d)+
  geom_point() +
  scale_x_log10()+
  scale_y_log10()+
  stat_smooth(method = lm)

library(smatr)
g <- sma(bvtv~mass_g, data = d, slope.test = 0, log = "xy", robust = T)

print(g) # slope = 0.315, p = 2.22e-16: positive allometry (slope is > 0). It's right around what the slope would be for isometry in a linear measurement, which is interesting, but there's no reason for it to be that because it is a dimensionless ratio. Bone Volume(mm^3)/Total volume (mm^3): the units cancel out. Right?

plot(g)
```

```{r}
ggplot(aes(x=tbth, y = mass_g), data = d)+
  geom_point() +
  scale_x_log10()+
  scale_y_log10()+
  stat_smooth(method = lm)

g2 <- sma(tbth~mass_g, data = d, slope.test = 1/3, log = "xy", robust = T)
print(g2) # slope = 0.328, p = 0.65: isometry (slope is not different from 1/3).  This differs from what Mielke et al 2018 uses because her body mass proxy is linear (she uses isometry = 1).

plot(g2)
```

```{r}

ggplot(aes(x=mass_g, y = tbsp), data = d)+
  geom_point() +
  scale_x_log10()+
  scale_y_log10()+
  geom_abline(slope = 0.271, intercept = -2.02)

g3 <- sma(tbsp~mass_g, data = d, slope.test = 1/3, log = "xy", robust = T)
print(g3) # slope = 0.271, p = 0.051... Isometry? Right on the damn edge. If it's negative allometry (slope is less than 1/3 is the test), it's very slight. This differs from what Mielke et al 2018 uses because her body mass proxy is linear (she uses isometry = 1).
plot(g3)
```

```{r}
ggplot(aes(x=connd, y = mass_g), data = d)+
  geom_point() +
  scale_x_log10()+
  scale_y_log10()+
  #stat_smooth(method = lm)+
  geom_abline(slope = -1, intercept = 5.931)

g4 <- sma(connd~mass_g, data = d, slope.test = -1, log = "xy", robust = T) # isometric slope -1 because it's connectivity per unit volume (1/mm^3). This differs from what Mielke et al 2018 uses because her body mass proxy is linear (she uses -3).
print(g4) #slope = -1, p = 0.99: Isometry (slope not diff from -1)
plot(g4)
```