---
title: "Context-dependent evolution of high trophic position drives functional disparity in subterranean crustaceans"
author: "Ester Premate, Špela Borko, Florian Altermatt, Cene Fišer"
date: "January 2023"
output:
  html_document:
    toc: TRUE
    toc_depth: 5
    toc_float: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, message=FALSE, warning=FALSE}
library(readxl)
library(phytools)
library(sensiPhy)
library(phylolm)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(car)
library(mvMORPH)
library(ellipse)
library(BAT)
library(ggtree)
library(stringr)
```

### 1) IMPORT & PREPARE DATA

##### Set working directory
```{r}
#set working directory
#setwd()
```

##### Load data
```{r}
#load Niphargus trait data
trait_data <- read_excel("Niphargus_morpho_data.xlsx", 
    col_types = c("text", "text", "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"), na = "NA")

head(trait_data)

#load Niphargus phylogenetic tree
tree <- read.nexus(file = "Niphargus_tree.nex")
```

##### Mutate and filter data
```{r}
#compute gnathopod perimeters from three gnathopod measures (sum of gnathopod length, width and diagonal)
trait_data <- trait_data %>% rowwise() %>% mutate(gnI_perim = sum(gnI_l, gnI_w, gnI_d))
trait_data <- trait_data %>% rowwise() %>% mutate(gnII_perim = sum(gnII_l, gnII_w, gnII_d))
#remove excess columns
trait_data <- trait_data[,-c(12:17)]
#remove species with missing morphological data
#N gallicus, N patrizii, N virei C
trait_data <- trait_data %>% drop_na(c(4:13))
```

### 2) DATA ANALYSIS 1: PHYLOGENETIC PCA

#### 2.1) Phylogenetic regression
Run phylogenetic regression: obtain residuals trait ~ body length
```{r}
#remove species with unknown ecology
trait_data <- trait_data %>% drop_na(c(2))
#prune tree
species <- as.vector(trait_data$species)
tree <- keep.tip(phy = tree, tip = species)
#assign rownames
trait_data <- data.frame(trait_data, row.names = 1)

#iterate phylogenetic regression function through columns 
phylolm <- lapply(trait_data[,-c(1:2)], function(y){
  
        phylolm(formula = y ~ trait_data$body_length, 
                data = trait_data, phy = tree, 
                model = "BM", track = FALSE)
  
})

#extract residuals from lists
residuals <- lapply(phylolm, function (x) x[c("residuals")])

#convert to dataframe and assign original column names
residuals_df <- as.data.frame(residuals)
colnames(residuals_df) <- names(residuals)
#add back scaled body length column
body_length <- scale(trait_data$body_length)
residuals_df <- cbind(residuals_df, body_length)

rm(residuals, phylolm, body_length)
```

#### 2.2) Phylogenetic PCA

Morphological traits were divided in three groups; "body shape" includes coxae II and III lengths and pereopod V-VII basis widths; "locomotion" includes pereopod V-VII lengths and body length, and "trophic position" includes gnathopod I and II perimeters.

```{r}
#create data matrices for phylogenetic PCA 
shape <- as.matrix(residuals_df[,c(4:8)]) #Body shape
locom <- as.matrix(residuals_df[,c(1:3, 11)]) #Locomotion
tp <- as.matrix(residuals_df[,c(9,10)])#Trophic position
```

```{r}
#run PCAs
shape_pca <- phyl.pca(tree = tree, shape, method = "BM", mode = "cov") 
locom_pca <- phyl.pca(tree = tree, locom, method = "BM", mode = "cov") 
tp_pca <- phyl.pca(tree = tree, tp, method = "BM", mode = "cov")

#save PC1 scores
pc1_scores <- as.data.frame(cbind(locom_pca$S[,1], shape_pca$S[,1], tp_pca$S[,1]))
colnames(pc1_scores) <- c("PC1_locom", "PC1_shape", "PC1_tp")
```

##### Body shape PCA results
```{r}
#PCA diagnostics
summary(shape_pca) #PC1: explains 71% of variance
shape_pca$L #Loadings: PC1: negatively correlates with all parameters.
```
##### Locomotion PCA results
```{r}
#PCA diagnostics
summary(locom_pca) #PC1: explains 82% of variance
locom_pca$L #Loadings: PC1: positively correlates with all parameters except body length.
```
##### Trophic position PCA results
```{r}
#PCA diagnostics
summary(tp_pca) #PC1: explains 94% of variance
tp_pca$L #Loadings: PC1: positively correlates with all parameters.
```
```{r}
rm(shape, shape_pca, locom, locom_pca, tp, tp_pca)
```


### 3) DATA ANALYSIS 2: REGRESSIONS BETWEEN PCs 
```{r}
#add ecology column to pc1_scores dataframe
ecology <- as.data.frame(trait_data[,1])
colnames(ecology) <- "ecology"
rownames(ecology) <- rownames(trait_data)
pc1_scores <- merge(pc1_scores, ecology, by = "row.names")
rownames(pc1_scores) <- pc1_scores[,1]
pc1_scores[,1] <- NULL
```

#### 3.1) Interstitial
##### Regression 1: PC1 tp ~ PC1 locomotion
```{r}
int_reg1 <- lm(PC1_tp ~ PC1_locom, data = pc1_scores %>% filter(ecology == "Interstitial"))
#plot diagnostics
par(mfrow = c(2, 2))
plot(int_reg1)
```

##### Regression 2: PC1 tp ~ PC1 shape
```{r}
int_reg2 <- lm(PC1_tp ~ PC1_shape, data = pc1_scores %>% filter(ecology == "Interstitial"))
#plot diagnostics
par(mfrow = c(2, 2))
plot(int_reg2) 
```

#### 3.2) Cave lake
##### Regression 1: PC1 tp ~ PC1 locomotion
```{r}
lake_reg1 <- lm(PC1_tp ~ PC1_locom, data = pc1_scores %>% filter(ecology == "Cave lake"))
#plot diagnostics
par(mfrow = c(2, 2))
plot(lake_reg1) 
```

##### Regression 2: PC1 tp ~ PC1 shape
```{r}
lake_reg2 <- lm(PC1_tp ~ PC1_shape, data = pc1_scores %>% filter(ecology == "Cave lake"))
#plot diagnostics
par(mfrow = c(2, 2))
plot(lake_reg2) 
```

#### 3.3) Cave stream
##### Regression 1: PC1 tp ~ PC1 locomotion
```{r}
stream_reg1 <- lm(PC1_tp ~ PC1_locom, data = pc1_scores %>% filter(ecology == "Cave stream"))
#plot diagnostics
par(mfrow = c(2, 2))
plot(stream_reg1) 
```

##### Regression 2: PC1 tp ~ PC1 shape
```{r}
stream_reg2 <- lm(PC1_tp ~ PC1_shape, data = pc1_scores %>% filter(ecology == "Cave stream"))
#plot diagnostics
par(mfrow = c(2, 2))
plot(stream_reg2) 
```

### 4) CORRELATIONS PLOTS & LM SUMMARIES

LM models summaries below each plot.

##### Interstitial species 
Plot 1
```{r}
theme_corr_plots <- theme_bw() + theme(axis.text = element_text(color = "black", size = 10),
          axis.title = element_text(color = "black", size = 10),
          panel.grid.major = element_blank(), panel.grid.minor = element_blank())
```

```{r}
ggplot(data = pc1_scores %>% filter(ecology == "Interstitial"), 
       mapping = aes(x = PC1_locom, y = PC1_tp)) +
    geom_point(color = "#ff7f00", size = 3, alpha = 0.8) +
    geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5) +
    ylab("PC1 Trophic apparatus") +
    xlab("PC1 Locomotion") +
    scale_y_continuous(labels = label_number(accuracy = 0.1)) +
    scale_x_continuous(labels = label_number(accuracy = 0.1)) +
    theme_corr_plots
#ggsave("inter1.png", width = 8, height = 8, units = "cm")
```

LM summary
```{r, echo=FALSE }
summary(int_reg1)
```

Plot2
```{r}
ggplot(data = pc1_scores %>% filter(ecology == "Interstitial"), 
       mapping = aes(x = PC1_shape, y = PC1_tp)) +
    geom_point(color = "#ff7f00", size = 3, alpha = 0.8) +
    geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5) +
    ylab("PC1 Trophic apparatus") +
    xlab("PC1 Body shape")  +
    scale_y_continuous(labels = label_number(accuracy = 0.1)) +
    scale_x_continuous(labels = label_number(accuracy = 0.1)) +
    theme_corr_plots
#ggsave("inter2.png", width = 8, height = 8, units = "cm")
```

LM summary
```{r, echo=FALSE}
summary(int_reg2)
```

##### Cave lake species 
Plot 1
```{r}
ggplot(data = pc1_scores %>% filter(ecology == "Cave lake"), 
       mapping = aes(x = PC1_locom, y = PC1_tp)) +
    geom_point(color = "#4dae48", size = 3, alpha = 0.8) +
    geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5) +
    ylab("PC1 Trophic apparatus") +
    xlab("PC1 Locomotion")  +
    scale_y_continuous(labels = label_number(accuracy = 0.1)) +
    scale_x_continuous(labels = label_number(accuracy = 0.1)) +
    theme_corr_plots
#ggsave("lake1.png", width = 8, height = 8, units = "cm")
```

LM summary
```{r,echo = FALSE}
summary(lake_reg1)
```

Plot 2
```{r}
ggplot(data = pc1_scores %>% filter(ecology == "Cave lake"), 
       mapping = aes(x = PC1_shape, y = PC1_tp)) +
    geom_point(color = "#4dae48", size = 3, alpha = 0.8) +
    geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5, linetype = "dashed") +
    ylab("PC1 Trophic apparatus") +
    xlab("PC1 Body shape")  +
    scale_y_continuous(labels = label_number(accuracy = 0.1)) +
    scale_x_continuous(labels = label_number(accuracy = 0.1)) +
    theme_corr_plots
#ggsave("lake2.png", width = 8, height = 8, units = "cm")
```

LM summary
```{r, echo=FALSE}
summary(lake_reg2)
```

##### Cave stream species 
Plot 1
```{r}
ggplot(data = pc1_scores %>% filter(ecology == "Cave stream"), 
       mapping = aes(x = PC1_locom, y = PC1_tp)) +
    geom_point(color = "#7570b2", size = 3, alpha = 0.8) +
    geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5) +
    ylab("PC1 Trophic apparatus") +
    xlab("PC1 Locomotion")  +
    scale_y_continuous(labels = label_number(accuracy = 0.1)) +
    scale_x_continuous(labels = label_number(accuracy = 0.1)) +
    theme_corr_plots
#ggsave("stream1.png", width = 8, height = 8, units = "cm")
```

LM summary
```{r, echo=FALSE}
summary(stream_reg1)
```

Plot 2
```{r}
ggplot(data = pc1_scores %>% filter(ecology == "Cave stream"), 
       mapping = aes(x = PC1_shape, y = PC1_tp)) +
    geom_point(color = "#7570b2", size = 3, alpha = 0.8) +
    geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5) +
    ylab("PC1 Trophic apparatus") +
    xlab("PC1 Body shape")  +
    scale_y_continuous(labels = label_number(accuracy = 0.1)) +
    scale_x_continuous(labels = label_number(accuracy = 0.1)) +
    theme_corr_plots #+
#ggsave("stream2.png", width = 8, height = 8, units = "cm")
```

LM summary
```{r,echo = FALSE}
summary(stream_reg2)
```
```{r}
#remove objects
rm(list = ls(pattern = "reg"))
```

### 5) DATA ANALYSIS 3: FUNCTIONAL DIVERSITY
#### 5.1) Functional diversity of real data and data simulated from phylogeny

##### Prepare tree and ecology data
```{r}
#tree
tree <- ladderize(tree, right = FALSE)

#ecological data (categorical)
ecology_ggtree <- cbind(rownames(ecology), data.frame(ecology, row.names=NULL))
colnames(ecology_ggtree) <- c("id","ecology")
ecology_ggtree <- ecology_ggtree[match(tree$tip.label, ecology_ggtree$id),]  

ecology_v <- setNames(ecology$ecology, rownames(ecology))
ecology_v <- ecology_v[tree$tip.label]
ecology_v <- as.factor(ecology_v)

# simmap tree
tree_simmap <- make.simmap(tree, ecology_v, model = "ER", nsim = 1)
plot(tree_simmap)
```

##### Run mvMORPH with different evolutionary models 
```{r}
#mvMorph 
model_BM1 <- mvBM(tree_simmap, pc1_scores[,1:3], model = "BM1")
model_BMM <- mvBM(tree_simmap,pc1_scores[,1:3], model = "BMM")
model_EB <- mvEB(tree,pc1_scores[,1:3])
model_OU1 <- mvOU(tree_simmap,pc1_scores[,1:3], model = "OU1")
model_OUM <- mvOU(tree_simmap,pc1_scores[,1:3], model = "OUM")
# merge into list
models <- list(model_BM1,model_BMM,model_EB,model_OU1,model_OUM)
#choose the best model
aic <- aicw(models)
aic
```

```{r}
#run simulations with the best model (OUM)
simul_phylo <- mvSIM(tree_simmap, nsim = 999, param = model_OUM)
```

```{r}
rm(list = ls(pattern = "model"))
rm(aic)
```

```{r}
#check first df
head(simul_phylo[[1]])
```

##### Construct functional trees 
First try test tree to see which method is the best (simulated data)
```{r}
testtree <- tree.build(trait = simul_phylo[[1]], distance = "euclidean", func = "best")
testtree
#the best tree is given by neighbor-joining algorithm
```
```{r}
ggtree(testtree, layout = "rectangular", branch.length = "none") + geom_tiplab(size = 1.5)
```

Construct 1000 functional trees from simulated data
```{r}
simul_func_trees <- lapply(simul_phylo, function(x){
  tree.build(trait = x, distance = "euclidean", func = "nj")
})
```

Construct functional tree from real data
```{r}
real_func_tree <- tree.build(trait = pc1_scores[,1:3], distance = "euclidean", func = "nj")
```

Create community x species matrix
```{r}
ecology$species <- rownames(ecology)
row.names(ecology) <- NULL
ecology <- ecology[,c(2,1)]
comm_sp_matrix <- reshape2::acast(ecology, ecology ~ species)
comm_sp_matrix[is.na(comm_sp_matrix)] <- 0
comm_sp_matrix[comm_sp_matrix == "Interstitial" | comm_sp_matrix == "Cave stream" | comm_sp_matrix == "Cave lake"] <- 1
#add another row to a matrix; "community" comprising all species
All <- rep(1,185)
#bind matrix and vector together
comm_sp_matrix <- rbind(comm_sp_matrix, All)
class(comm_sp_matrix) <- "numeric"
rm(All)
```

##### Alpha diversity
Calculate simulated alpha functional diversity across trees
```{r}
#check if matrix column names and tip labels are the same
setdiff(colnames(comm_sp_matrix),simul_func_trees[[1]]$tip.label)
```

```{r}
#calculate alpha diversity
alpha_sim1 <- lapply(simul_func_trees, function(x) {alpha(comm = comm_sp_matrix, tree = x)})
```

```{r}
#check first element
alpha_sim1[[1]]
```

```{r}
#extract richness from list for every tree
alpha_sim_df1 <- as.data.frame(do.call(rbind, alpha_sim1))
alpha_sim_df1 <- tibble::rownames_to_column(alpha_sim_df1, "habitat")
```


```{r}
#replace values
alpha_sim_df1$habitat <- str_replace_all(alpha_sim_df1$habitat,
                                         c("Cave.lake.\\d+" = "Cave lake",
                                           "Cave.lake" = "Cave lake",
                                           "Interstitial.\\d+" = "Interstitial",
                                           "Cave.stream.\\d+" = "Cave stream",
                                           "Cave.stream" = "Cave stream",
                                           "All.\\d+" = "All"))
```

```{r}
unique(alpha_sim_df1$habitat)
```


Calculate real alpha functional diversity 
```{r}
#calculate alpha diversity
alpha_real1 <- as.data.frame(alpha(comm = comm_sp_matrix, tree = real_func_tree))
alpha_real1 <- tibble::rownames_to_column(alpha_real1, "habitat")
```

Calculate p values
```{r}
#p-value cave lake phylo null model
ses(alpha_real1[1,2], (alpha_sim_df1 %>% filter(habitat == "Cave lake"))$Richness)

#p-value cave stream phylo null model
ses(alpha_real1[2,2], (alpha_sim_df1 %>% filter(habitat == "Cave stream"))$Richness)

#p-value interstitial phylo null model
ses(alpha_real1[3,2], (alpha_sim_df1 %>% filter(habitat == "Interstitial"))$Richness)

#p-value all phylo null model
ses(alpha_real1[4,2], (alpha_sim_df1 %>% filter(habitat == "All"))$Richness)
```

##### Beta diversity 
```{r}
#drop row "all" from comm_sp_matrix 
comm_sp_matrix_beta <- comm_sp_matrix[-4,]
```

```{r}
#calculate beta diversity
beta_sim1 <- lapply(simul_func_trees, function(x) {beta(comm = comm_sp_matrix_beta, tree = x, func = "jaccard", abund = FALSE)})
```

```{r}
#check first element
beta_sim1[[1]]
```

```{r}
beta_sim1[[1]]$Btotal
```

Extract data from lists
Btotal1: pair cave stream-cave lake
Btotal2: pair interstitial-cave lake
Btotal3: pair interstitial-cave stream
```{r}
#unlist
beta_sim_df1 <- lapply(beta_sim1, function(x){unlist(x)})
#get data for each pair of habitats
#cave stream-interstitial
stream_interstitial <- as.data.frame(sapply(beta_sim_df1, "[[", "Btotal3"))
colnames(stream_interstitial) <- "beta_total"
stream_interstitial$comparison <- "Cave stream - Interstitial"
#cave stream-lake
stream_lake <- as.data.frame(sapply(beta_sim_df1, "[[", "Btotal1"))
colnames(stream_lake) <- "beta_total"
stream_lake$comparison <- "Cave stream - Cave lake"
#interstitial-lake
interstitial_lake <- as.data.frame(sapply(beta_sim_df1, "[[", "Btotal2"))
colnames(interstitial_lake) <- "beta_total"
interstitial_lake$comparison <- "Interstitial - Cave lake"
#bind together
beta_sim_df1 <- rbind(interstitial_lake, stream_lake, stream_interstitial)
```

Average simulated beta diversity across communities
```{r}
beta_sim_multi1 <- lapply(simul_func_trees, function(x){beta.multi(comm = comm_sp_matrix_beta, tree = x, func = "jaccard", abund = FALSE)})
```

```{r}
#check first element
beta_sim_multi1[[1]]
```

```{r}
#extract data
beta_sim_multi_df1 <- as.data.frame(do.call(rbind, beta_sim_multi1))
beta_sim_multi_df1 <- tibble::rownames_to_column(beta_sim_multi_df1, "beta_total")
beta_sim_multi_df1 <- beta_sim_multi_df1[,c(1,2)] %>%  filter(str_detect(beta_total, "Btotal"))
beta_sim_multi_df1$comparison <- "Average (all)"
beta_sim_multi_df1$beta_total <- NULL
colnames(beta_sim_multi_df1) <- c("beta_total", "comparison")
```

```{r}
#bind together with other pairs
beta_sim_df1 <- rbind(beta_sim_df1, beta_sim_multi_df1)
rm(beta_sim_multi_df1, interstitial_lake, stream_lake, stream_interstitial)
```

Calculate real beta functional diversity 
```{r}
#calculate beta diversity
beta_real1 <- beta(comm = comm_sp_matrix_beta, tree = real_func_tree, func = "jaccard", abund = FALSE)
beta_real1
```

Average real beta diversity across communities
```{r}
beta_real_multi1 <- beta.multi(comm = comm_sp_matrix_beta, tree = real_func_tree, func = "jaccard", abund = FALSE)
```

```{r}
#keep only total and average
beta_real_multi1 <- as.data.frame(beta_real_multi1[1,1])
colnames(beta_real_multi1) <- "beta_total"
beta_real_multi1$comparison <- "Average (all)"
```

```{r}
beta_real1 <- as.data.frame(unlist(beta_real1))
beta_real1 <- tibble::rownames_to_column(beta_real1, "beta_total")[1:3,]
beta_real1$comparison <- c("Cave stream - Cave lake", "Interstitial - Cave lake", "Cave stream - Interstitial")
beta_real1$beta_total <- NULL
colnames(beta_real1) <- c("beta_total", "comparison")
```

```{r}
#bind together with other pairs
beta_real_df1 <- rbind(beta_real1, beta_real_multi1)
beta_real_df1$beta_total <- as.numeric(beta_real_df1$beta_total)
```

Calculate p values
```{r}
#p-value stream-lake phylo null model
ses(beta_real_df1[1,1], (beta_sim_df1 %>% filter(comparison == "Cave stream - Cave lake"))$beta_total)

#p-value interstitial-lake phylo null model
ses(beta_real_df1[2,1], (beta_sim_df1 %>% filter(comparison == "Interstitial - Cave lake"))$beta_total)

#p-value stream-interstitial phylo null model
ses(beta_real_df1[3,1], (beta_sim_df1 %>% filter(comparison == "Cave stream - Interstitial"))$beta_total)

#p-value all phylo null model
ses(beta_real_df1[4,1], (beta_sim_df1 %>% filter(comparison == "Average (all)"))$beta_total)
```

```{r}
rm(alpha_sim1, beta_real_multi1, beta_sim1, beta_sim_multi1, real_func_tree, simul_func_trees, simul_phylo, tree_simmap, ecology_ggtree, testtree)
```

#### 5.2) Functional diversity of real data and data randomly drawn from real data

##### Create simulations from real data
```{r}
#edit pc1_scores df
pc1_scores$species <- rownames(pc1_scores)
row.names(pc1_scores) <- NULL
pc1_scores <- pc1_scores[,c(5, 1:4)]
#create empty list
simul_list <- vector("list", length = 999)
#fill with the same dataframe
simul_list <- lapply(simul_list, function(x) {pc1_scores})
```

```{r}
#randomize dataframes (keep fixed ecological groups)
simul_list <- lapply(simul_list, function(x){
  x <- x %>% group_by(ecology) %>% mutate(locom_sim = sample(PC1_locom)) %>% 
    mutate(shape_sim = sample(PC1_shape)) %>% 
    mutate(tp_sim = sample(PC1_tp))
})
#convert to dataframes
simul_list <- lapply(simul_list, function(x) {as.data.frame(x)})
```

```{r}
#check
head(simul_list[[1]])
```

```{r}
#species as rownames
simul_list <- lapply(simul_list, function(x) {rownames(x) <- x$species;x})
#remove unneccessary columns: species, real values and ecology
simul_list <- lapply(simul_list, function(x) {x <- x[,c(6:8)]})
#check
head(simul_list[[1]])
```

##### Construct functional trees
First try test tree to see which method is best
```{r}
testtree <- tree.build(trait = simul_list[[1]], distance = "euclidean", func = "best")
testtree #the best tree is given by neighbor-joining algorithm
```
Construct functional tree from simulated data
```{r}
simul_func_trees <- lapply(simul_list, function(x){
  tree.build(trait = x, distance = "euclidean", func = "nj")
})
```

Construct functional tree from real data
```{r}
rownames(pc1_scores) <- pc1_scores$species
pc1_scores <- pc1_scores[,c(-1,-5)]
real_func_tree <- tree.build(trait = pc1_scores, distance = "euclidean", func = "nj")
```

##### Alpha diversity
Calculate simulated alpha functional diversity across trees
```{r}
#check if matrix column names and tip labels are the same
setdiff(colnames(comm_sp_matrix),simul_func_trees[[1]]$tip.label)
```

```{r}
#calculate alpha diversity
alpha_sim2 <- lapply(simul_func_trees, function(x) {alpha(comm = comm_sp_matrix, tree = x)})
```

```{r}
#check first element
alpha_sim2[[1]]
```

```{r}
#extract richness from list for every tree
alpha_sim_df2 <- as.data.frame(do.call(rbind, alpha_sim2))
alpha_sim_df2 <- tibble::rownames_to_column(alpha_sim_df2, "habitat")
```


```{r}
#replace values
alpha_sim_df2$habitat <- str_replace_all(alpha_sim_df1$habitat,
                                         c("Cave.lake.\\d+" = "Cave lake",
                                           "Cave.lake" = "Cave lake",
                                           "Interstitial.\\d+" = "Interstitial",
                                           "Cave.stream.\\d+" = "Cave stream",
                                           "Cave.stream" = "Cave stream",
                                           "All.\\d+" = "All"))
```


```{r}
unique(alpha_sim_df2$habitat)
```

Calculate p values
```{r}
#p-value cave lake phylo null model
ses(alpha_real1[1,2], (alpha_sim_df2 %>% filter(habitat == "Cave lake"))$Richness)

#p-value cave stream phylo null model
ses(alpha_real1[2,2], (alpha_sim_df2 %>% filter(habitat == "Cave stream"))$Richness)

#p-value interstitial phylo null model
ses(alpha_real1[3,2], (alpha_sim_df2 %>% filter(habitat == "Interstitial"))$Richness)

#p-value all phylo null model
ses(alpha_real1[4,2], (alpha_sim_df2 %>% filter(habitat == "All"))$Richness)
```
##### Beta diversity
Calculate simulated beta functional diversity 
```{r}
#drop row "all" from comm_sp_matrix 
comm_sp_matrix_beta <- comm_sp_matrix[-4,]
```

```{r}
#calculate beta diversity
beta_sim2 <- lapply(simul_func_trees, function(x) {beta(comm = comm_sp_matrix_beta, tree = x, func = "jaccard", abund = FALSE)})
```

```{r}
#check first element
beta_sim2[[1]]
```

Extract data from lists
Btotal1: pair cave stream-cave lake
Btotal2: pair interstitial-cave lake
Btotal3: pair interstitial-cave stream
```{r}
#unlist
beta_sim_df2 <- lapply(beta_sim2, function(x){unlist(x)})
#get data for each pair of habitats
#cave stream-interstitial
stream_interstitial <- as.data.frame(sapply(beta_sim_df2, "[[", "Btotal3"))
colnames(stream_interstitial) <- "beta_total"
stream_interstitial$comparison <- "Cave stream - Interstitial"
#cave stream-lake
stream_lake <- as.data.frame(sapply(beta_sim_df2, "[[", "Btotal1"))
colnames(stream_lake) <- "beta_total"
stream_lake$comparison <- "Cave stream - Cave lake"
#interstitial-lake
interstitial_lake <- as.data.frame(sapply(beta_sim_df2, "[[", "Btotal2"))
colnames(interstitial_lake) <- "beta_total"
interstitial_lake$comparison <- "Interstitial - Cave lake"
#bind together
beta_sim_df2 <- rbind(stream_interstitial, stream_lake, interstitial_lake)
```

Average simulated beta diversity across communities
```{r}
beta_sim_multi2 <- lapply(simul_func_trees, function(x){beta.multi(comm = comm_sp_matrix_beta, tree = x, func = "jaccard", abund = FALSE)})
```

```{r}
#check first element
beta_sim_multi2[[1]]
```

```{r}
#extract data
beta_sim_multi_df2 <- as.data.frame(do.call(rbind, beta_sim_multi2))
beta_sim_multi_df2 <- tibble::rownames_to_column(beta_sim_multi_df2, "beta_total")
beta_sim_multi_df2 <- beta_sim_multi_df2[,c(1,2)] %>%  filter(str_detect(beta_total, "Btotal"))
beta_sim_multi_df2$comparison <- "Average (all)"
beta_sim_multi_df2$beta_total <- NULL
colnames(beta_sim_multi_df2) <- c("beta_total", "comparison")
```

```{r}
#bind together with other pairs
beta_sim_df2 <- rbind(beta_sim_df2, beta_sim_multi_df2)
rm(beta_sim_multi_df2, stream_interstitial, stream_lake, interstitial_lake)
```

Calculate p values
```{r}
#p-value stream-lake phylo null model
ses(beta_real_df1[1,1], (beta_sim_df2 %>% filter(comparison == "Cave stream - Cave lake"))$beta_total)

#p-value interstitial-lake phylo null model
ses(beta_real_df1[2,1], (beta_sim_df2 %>% filter(comparison == "Interstitial - Cave lake"))$beta_total)

#p-value stream-interstitial phylo null model
ses(beta_real_df1[3,1], (beta_sim_df2 %>% filter(comparison == "Cave stream - Interstitial"))$beta_total)

#p-value all phylo null model
ses(beta_real_df1[4,1], (beta_sim_df2 %>% filter(comparison == "Average (all)"))$beta_total)
```




#### Plots 

Alpha diversity
```{r}
#combine data in single df
alpha_sim_df1$simul <- "Phylogenetic null model"
alpha_sim_df2$simul <- "Non-phylogenetic null model"
alpha_sim_df <- rbind(alpha_sim_df1, alpha_sim_df2)
```


```{r}
quantiles_a_phy <- alpha_sim_df1 %>% group_by(habitat) %>% 
  summarise(int = quantile(Richness, c(0.025, 0.975)))

quantiles_a_nonphy <- alpha_sim_df2 %>% group_by(habitat) %>% 
  summarise(int = quantile(Richness, c(0.025, 0.975))) 
```

```{r}
ggplot() + 
  geom_density(data = alpha_sim_df %>% filter(simul == "Phylogenetic null model"), 
              mapping = aes(x = Richness, fill = habitat), 
              alpha = 0.5, color = NA) +
  geom_density(data = alpha_sim_df %>% filter(simul == "Non-phylogenetic null model"), 
               mapping = aes(x = Richness, color = habitat), alpha = 1) +
  geom_point(data = alpha_real1, 
             mapping = aes(x = Richness, y = 0, fill = habitat), shape = 21, size = 3) + 

  scale_fill_manual(values = c("grey40", "#4dad47", "#7570b1","#f37d20")) +
  scale_color_manual(values = c("grey40", "#4dad47", "#7570b1","#f37d20")) +

  facet_wrap(~habitat, scales = "fixed", nrow = 4, ncol = 1) +
  geom_vline(data = quantiles_a_phy, aes(xintercept = int, color = habitat), alpha = 0.7) +
  geom_vline(data = quantiles_a_nonphy, aes(xintercept = int, color = habitat), 
             linetype = "dashed") +
  
  xlab(expression(alpha*-functional~diversity~(richness))) +
  ylab("Density") +
  theme_bw() +
  theme(axis.text = element_text(color = "black"),
        legend.position = "none",
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill = "white"))
  
#ggsave("sim_alpha_density.svg", plot = last_plot(), dpi = 300, width = 4, height = 8)
```

Beta diversity
```{r}
#combine data in single df
beta_sim_df1$simul <- "Phylogenetic null model"
beta_sim_df2$simul <- "Non-phylogenetic null model"
beta_sim_df <- rbind(beta_sim_df1, beta_sim_df2)
```

```{r}
quantiles_b_phy <- beta_sim_df1 %>% group_by(comparison) %>% 
  summarise(int = quantile(beta_total, c(0.025, 0.975)))
quantiles_b_nonphy <- beta_sim_df2 %>% group_by(comparison) %>% 
  summarise(int = quantile(beta_total, c(0.025, 0.975))) 
```

```{r}
ggplot() + 
  geom_density(data = beta_sim_df %>% filter(simul == "Phylogenetic null model"), 
              mapping = aes(x = beta_total), alpha = 0.5, fill = "grey40", color = NA) +
  
  geom_density(data = beta_sim_df %>% filter(simul == "Non-phylogenetic null model"), 
               mapping = aes(x = beta_total), alpha = 1) +
  geom_point(data = beta_real_df1, 
             mapping = aes(x = beta_total, y = 0), shape = 21, size = 3) + 

 facet_wrap(~comparison, scales = "fixed", nrow = 4, ncol = 1) +
  
  geom_vline(data = quantiles_b_phy, aes(xintercept = int), alpha = 0.7, color = "grey40") +
  geom_vline(data = quantiles_b_nonphy, aes(xintercept = int), linetype = "dashed", color = "black") +
  
  xlab(expression(beta*-functional~diversity~(total))) +
  ylab("Density") +
  theme_bw() +
  theme(axis.text = element_text(color = "black"),
        legend.position = "none",
        panel.grid.major.y = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill = "white"))
 
#ggsave("sim_beta_density.svg", plot = last_plot(), dpi = 300, width = 4, height = 8)
#colors to this plot were added in Adobe Illustrator
```

```{r, echo = FALSE}
#package versions
sessionInfo()
```



```{r, eval=FALSE, echo = FALSE, results= "hide"}
fig1 <- read_excel("Fig1,2/fig1data.xlsx")
```

```{r, eval=FALSE, echo = FALSE, results= "hide"}
ggplot() + geom_point(data = fig1, mapping = aes(x = Trait1, y = Trait2)) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        axis.text = element_blank(),
        panel.border = element_blank(),
        axis.line = element_line(arrow = arrow(length = unit(10,'pt'))),
        axis.ticks = element_blank(),
        axis.title = element_text(size = 18)) +
  xlab("Trait 1") +
  ylab("Trait 2")
#ggsave("fig1_base.png", plot = last_plot(), dpi = 300, width = 4, height = 4)
```

