---
title: "Cichlid v Centrarchid Integration_EVO-23-0246"
author: "A. Roberts-Hugghis with input from E.Burress"
date: "2023-06-27"
output: html_document
---

## R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>. When you click the **Knit** button, a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. 

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

# <--------------- SET WORKING DIRECTORY ---------------> #
setwd("~/Documents/CichvCent_Inter/")


# <--------------- LOAD PACKAGES ---------------> #
require(tidyverse) 
require(purrr)        
require(phytools) 
require(geiger)
require(geomorph)
require(caper)
require(ggplot2)
require(ggrepel)
require(rgl)
require(mvnormtest)
require(effsize)
require(varhandle)
require(motmot)
require(diversitree)
require(devtools)
require(grid)
require(gridExtra)
require(gtable)
require(reshape)
require(lattice)
require(cowplot)
#install_github("mlcollyer/RRPP")
require(RRPP) 
require(corrplot)   
require(hypervolume)    
require(deeptime)  
require(ggphylomorpho)
require(cowplot)


# <--------------- LOAD AND FORMAT DATA ---------------> #
#Load data frame + retrieve key information 
SR.data.mean <- read.csv("SR_data_mean.csv")

ComboFin2_new <- read.csv("ComboFin2_new.csv")
ComboFin2_new <- ComboFin2_new[,-1]
row.names(ComboFin2_new) <- ComboFin2_new$Species3
Species3 <- ComboFin2_new$Species3

#Load tree 
ComboTreeFin2 <- read.tree("ComboTreeFin2.tre")
plot(ComboTreeFin2)
```        

        
```{r PRINCIPAL COMPONENT ANALYSES, warning=FALSE}    

####################################
###### ALL SPECIES + TRAITS ########
####################################

# <--------------- FORMAT DATA ---------------> #
#Create data frame for PCA of all traits and all species; exclude columns with species and family info
prepca <- ComboFin2_new[,3:22]

#Make family a factor as it will be used as categorical variable. Assign colors to each species based on categorical variable.
Fam2 <- as.factor(ComboFin2_new$Fam2) 
clr2 <- c("Cichlidae" = "rosybrown", "Centrarchidae" = "skyblue4")
Fam2 <- fct_relevel(Fam2, c("Cichlidae","Centrarchidae"))


# <--------------- RUN + VISUALIZE PCA ---------------> #
#Perform the PCA using the prcomp function. Here, we specify that we are using the correlation matrix using "scale = TRUE" because the data has different units/scales.
combo.pca <- prcomp(prepca, scale =TRUE)

#Summarize PCA results 
summary(combo.pca) #Informs us of which axes explain the most variation in the data
plot(combo.pca, main = "Variance Explained by PCs", xlab = "Principal Components") #Creates a histogram of variance explained by PCs
biplot(combo.pca, cex = .5) #Shows species in morphospace and how variables load on which PCs. Can use choices to pick and choose variables to show you (e.g., biplot(combo.pca2.1, choices=c(4,5)))

#Creates a data frame of pc loadings (how much each variable is explained by each PC axis) and PC scores (location of each point (species) on each PC axis)
pc.loadings <- as.data.frame(combo.pca$rotation) 
pc.scores <- as.data.frame(combo.pca$x)

#Creates a data frame of two pc axes and the categorical variable you want to color by
gg1 <- data.frame('Fam' = ComboFin2_new$Fam, combo.pca$x[,1:2])
gg1$PC1 <- -(gg1$PC1) #Rotates axes
gg1$PC2 <- -(gg1$PC2) 

#Visualize PCA in ggplot
gmin1 <- ggplot(data = gg1, aes(x = PC1, y = PC2, col = Fam)) + 
  geom_point(size = 5, alpha=.85) +
  scale_colour_manual(values = clr2) +
  theme_classic() + 
  theme(legend.position="none", axis.text = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = .825) #+
  #geom_text(label=(ComboFin2_new$Species), size=1.5, hjust=0, vjust=0)
par(mar=c(0,0,1,0))
gmin1




#######################################################
###### CICHS & CENTS SEPARATELY - FIG 2 in PUB ########
#######################################################

# <--------------- FORMAT DATA ---------------> #
#Create data frame for PCA of cichs and cents separately 
ComboFin2_cich <- prepca[1:85, ]
ComboFin2_cent <- prepca[86:115, ]


# <--------------- RUN + VISUALIZE PCA ---------------> #
#Perform the PCA using the prcomp function. Here, we specify that we are using the correlation matrix using "scale = TRUE" because the data has different units/scales.
combo.pca2_cich <-prcomp(ComboFin2_cich, scale=TRUE) 
combo.pca2_cent <-prcomp(ComboFin2_cent, scale=TRUE) 


#Summarize PCA results 
summary(combo.pca2_cich) #Informs us of which axes explain the most variation in the data
summary(combo.pca2_cent) 

plot(combo.pca2_cich, main = "Variance Explained by PCs", xlab = "Principal Components") #Creates a histogram of variance explained by PCs
plot(combo.pca2_cent, main = "Variance Explained by PCs", xlab = "Principal Components")

biplot(combo.pca2_cich, cex = .5) #Shows species in morphospace and how variables load on which PCs. Can use choices to pick and choose variables to show you (e.g., biplot(combo.pca2.1, choices=c(4,5)))
biplot(combo.pca2_cent, cex = .5)


#Creates a data frame of pc loadings (how much each variable is explained by each PC axis) and PC scores (location of each point (species) on each PC axis)
pc.loadings_cich <- as.data.frame(combo.pca2_cich$rotation)
pc.loadings_cent <- as.data.frame(combo.pca2_cent$rotation)

pc.scores_cich <- as.data.frame(combo.pca2_cich$x)
pc.scores_cent <- as.data.frame(combo.pca2_cent$x)


#Creates a data frame of two pc axes and the categorical variable you want to color by
gg1_cich <- data.frame('Fam' = ComboFin2_new$Fam2[1:85], combo.pca2_cich$x[,1:2])
gg1_cich$PC1 <- -(gg1_cich$PC1)

gg1_cent <- data.frame('Fam' = ComboFin2_new$Fam2[86:115], combo.pca2_cent$x[,1:2])
gg1_cent$PC1 <- -(gg1_cent$PC1)


#Visualize PCA in ggplot
gmin1_cich <- ggplot(data = gg1_cich, aes(x = PC1, y = PC2, col = Fam)) + 
  geom_point(size = 5, alpha=.85) +
  scale_colour_manual(values = clr2) +
  theme_classic() + 
  theme(legend.position="none", axis.text = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = .825) #+
#geom_text(label=(ComboFin2$Species), size=1.5, hjust=0, vjust=0)
gmin1_cich 


gmin1_cent <- ggplot(data = gg1_cent, aes(x = PC1, y = PC2, col = Fam)) + 
  geom_point(size = 5, alpha=.85) +
  scale_colour_manual(values = clr2) +
  theme_classic() + 
  theme(legend.position="none", axis.text = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = .825) #+
#geom_text(label=(ComboFin2$Species), size=1.5, hjust=0, vjust=0)
gmin1_cent 




###################################################################
###### ORAL AND PHARYNGEAL JAWS SEPARATELY  - FIG 3 in PUB ########
###################################################################

# <--------------- FORMAT DATA ---------------> #
ComboFin2_OJ <- prepca[ ,1:13]
ComboFin2_PJ <-  prepca[ ,14:20]


# <--------------- RUN + VISUALIZE PCA ---------------> #
#Perform the PCA using the prcomp function. Here, we specify that we are using the correlation matrix using "scale = TRUE" because the data has different units/scales.
combo.pca2_OJ <-prcomp(ComboFin2_OJ, scale=TRUE) 
combo.pca2_PJ <-prcomp(ComboFin2_PJ, scale=TRUE) 


#Summarize PCA results 
summary(combo.pca2_OJ) #Informs us of which axes explain the most variation in the data
summary(combo.pca2_PJ)

plot(combo.pca2_OJ, main = "Variance Explained by PCs", xlab = "Principal Components") #Creates a histogram of variance explained by PCs
plot(combo.pca2_PJ, main = "Variance Explained by PCs", xlab = "Principal Components")

biplot(combo.pca2_OJ, cex = .5) #Shows species in morphospace and how variables load on which PCs. Can use choices to pick and choose variables to show you (e.g., biplot(combo.pca2.1, choices=c(4,5)))
biplot(combo.pca2_PJ, cex = .5)


#Creates a data frame of pc loadings (how much each variable is explained by each PC axis) and PC scores (location of each point (species) on each PC axis)
pc.loadings_OJ <- as.data.frame(combo.pca2_OJ$rotation)
pc.loadings_PJ <- as.data.frame(combo.pca2_PJ$rotation)

pc.scores_OJ <- as.data.frame(combo.pca2_OJ$x)
pc.scores_PJ <- as.data.frame(combo.pca2_PJ$x)


#Creates a data frame of two pc axes and the categorical variable you want to color by
gg1_OJ <- data.frame('Fam' = ComboFin2_new$Fam, combo.pca2_OJ$x[,1:2])
gg1_OJ$PC1 <- -(gg1_OJ$PC1)

gg1_PJ <- data.frame('Fam' = ComboFin2_new$Fam, combo.pca2_PJ$x[,1:2])


#Visualize PCA in ggplot
gmin1_OJ <- ggplot(data = gg1_OJ, aes(x = PC1, y = PC2, col = Fam)) + 
  geom_point(size = 5.25, alpha=.85) +
  scale_colour_manual(values = clr2) +
  theme_classic() + 
  theme(legend.position="none", axis.text = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = .825) #+
  #geom_text(label=(ComboFin2_OJ$Species), size=1.5, hjust=0, vjust=0)
gmin1_OJ 

gmin1_PJ <- ggplot(data = gg1_PJ, aes(x = PC1, y = PC2, col = Fam)) + 
  geom_point(size = 5, alpha=.85) +
  scale_colour_manual(values = clr2) +
  theme_classic() + 
  theme(legend.position="none", axis.text = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = .825) #+
   #geom_text(label=(ComboFin2_PJ$Species), size=1.5, hjust=0, vjust=0)
gmin1_PJ 
```




```{r ANALYSES OF VARIANCE, warning=FALSE}   

# <--------------- FORMAT DATA ---------------> #
#Create a matrix that has all the traits you want to analyze
integration_ANOVA_mat2 <- prepca %>%
  mutate_at(vars(Premax_Len3:UPJ_Dwteeth3), funs(scale(.)))
integration_ANOVA_mat3 <- as.matrix(integration_ANOVA_mat2)
OJintegration_ANOVA_mat <- integration_ANOVA_mat3 [,1:13]
PJintegration_ANOVA_mat <- integration_ANOVA_mat3 [,14:20]

#Create the x and y values for the ANOVA where x is the grouping variable and y is the matrix of traits we are evaluating
ANOVA_x2 <- ComboFin2_new$Fam2
names(ANOVA_x2) <- ComboFin2_new$Species3


# <--------------- RUN PHYLO MANOVAS ---------------> #
#All traits
inter_manova3 <- procD.pgls(integration_ANOVA_mat3~ANOVA_x2, phy = ComboTreeFin2, iter=9999)
summary(inter_manova3)

#Oral jaw traits
OJinter_manova <- procD.pgls(OJintegration_ANOVA_mat~ANOVA_x2, phy = ComboTreeFin2, iter=9999)
summary(OJinter_manova)

#Pharyngeal jaw traits
PJinter_manova <- procD.pgls(PJintegration_ANOVA_mat~ANOVA_x2, phy = ComboTreeFin2, iter=9999)
summary(PJinter_manova)


# <--------------- RUN PHYLO ANOVAS FOR EACH TRAIT ---------------> #
#Create an empty list for the results to go into
ANOVAresults  <- list()

#Run ANOVA for each trait using a for loop
for (i in 1:ncol(integration_ANOVA_mat2)) {
  ANOVAresults[i] <- procD.pgls(integration_ANOVA_mat3[,i]~ANOVA_x2, phy = ComboTreeFin2, iter=9999)
}
View(ANOVAresults)


# <--------------- RUN NON-PHYLO MANOVAS ---------------> #
#Can run non-phylo MANOVAs/ANOVAs to see if it’s the phylogeny that’s driving the result. If this is the case then you will see non-significant differences in mean trait values when running phylo MANOVAs/ANOVAs, but significant differences when running non-phylo MANOVAs/ANOVAs. This may be especially important when you see that your groups are in distinct spaces (i.e., showing little overlap) in a PCA, but not showing significant differences in mean trait values when running phylo MANOVAs/ANOVAs. If it IS  the phylogeny then it’s impossible to disentangle if the morphology, ecology, or other things are driving the morphological difference in the two groups.
inter_NONPHYmanova <- manova(integration_ANOVA_mat3~ANOVA_x2)
summary(inter_NONPHYmanova)
 

```




```{r DISPARITY, warning=FALSE}   

# <--------------- FORMAT DATA ---------------> #
disp_Combo3 <- as.matrix(prepca)
disp_OJ3 <- as.matrix(ComboFin2_OJ)
disp_PJ3 <- as.matrix(ComboFin2_PJ)

dispFAM <- as.vector(ComboFin2_new$Fam2)
names(dispFAM) <- ComboFin2_new$Species3


# <--------------- RUN QUICK MULTIVARIATE DISPARITY ANALYSES---------------> #
#Quick disparity analyses were run with avg.sq. This is average squared Euclidean distance among all pairs of points; this is the most common distance metric for disparity in macroevolution. This analysis provides no comparison of the groups!
(All_disarity3 <- disparity(phy=NULL, data = disp_Combo3, index = c("avg.sq")))
(OJ_disarity3 <- disparity(phy=NULL, data = disp_OJ3, index = c("avg.sq")))
(PJ_disarity3 <- disparity(phy=NULL, data = disp_PJ3, index = c("avg.sq"))) 


# <--------------- RUN MULTIVARIATE DISPARITY ANALYSES in GEOMORPH---------------> #
#Comparison between groups for all traits, OJs, and PJs in geomorph. This allows us to calculate the disparity in each of the two groups and compare it to determine if there are significant differences. No phylo included in these anaylses. 
#Results reported in Table 1
(morphoDF <- morphol.disparity(f1 = disp_Combo3 ~ dispFAM, groups = ~ dispFAM, iter = 9999, seed = NULL, print.progress = TRUE))

(morphoDF_OJ <- morphol.disparity(f1 = disp_OJ3 ~ dispFAM, groups = ~ dispFAM, iter = 9999, seed = NULL, print.progress = TRUE))

(morphoDF_PJ <- morphol.disparity(f1 = disp_PJ3 ~ dispFAM, groups = ~ dispFAM, iter = 9999, seed = NULL, print.progress = TRUE))


# <--------------- RUN UNIVARIATE DISPARITY ANALYSES in GEOMORPH---------------> #
#Set up data frame to catch disparity results
groups <- levels(as.factor(dispFAM))
trait_vec <- colnames(disp_Combo3)

TraitDisparity <- data.frame(matrix(ncol = ncol(disp_Combo3), nrow = length(groups)+2))
colnames(TraitDisparity) <- trait_vec
row.names(TraitDisparity) <- c(groups,"ratio","p-value")

#sequential_disparity
for(i in 1:ncol(disp_Combo3)){
  print(paste("Now running disparity for", trait_vec[i]))
  run_disparity <- morphol.disparity(f1 = disp_Combo3[,i] ~ dispFAM, groups = ~ dispFAM, partial = FALSE, 
                                     iter = 9999, seed = NULL, print.progress = F)
  
  #get pairwise disparities
  for (j in 1:length(groups)){
    TraitDisparity[j,i] <- run_disparity$Procrustes.var[j]
  }
  
  #get difference (multiples more or less)
  for (k in 1:ncol(TraitDisparity)){
    group1 <- TraitDisparity[1,k] #non-PHGY disparity
    group2 <- TraitDisparity[2,k] #PHGY disparity
    group_diff <- group1/group2
    TraitDisparity[3,k] <- group_diff
  }
  
  #get pairwise p-values on disparities
    TraitDisparity[4,i] <- run_disparity$PV.dist.Pval[1,2]
  
}


# <--------------- VISUALIZE DISPARITY RESULTS---------------> #

##########################################
###### BAR PLOT - NOT USED IN PUB ########
##########################################
#Retrieve variance info + format in df
morpholDisp <- as.data.frame(c(morphoDF$Procrustes.var,morphoDF_OJ$Procrustes.var,morphoDF_PJ$Procrustes.var))
names(morpholDisp) [1]  <- "Disparity"
morpholDisp$Group <- as.factor(c("Centrarchidae", "Cichlidae", "Centrarchidae", "Cichlidae", "Centrarchidae", "Cichlidae"))
morpholDisp$DispType<- as.factor(c("Overall", "Overall", "Oral Jaws", "Oral Jaws", "Pharyngeal Jaws", "Pharyngeal Jaws"))
morpholDisp$DispType <- fct_relevel(morpholDisp$DispType, c("Pharyngeal Jaws", "Oral Jaws", "Overall"))

#Compute std dev and add it to the df
morpholDisp$sd[1:2] <- sd(morpholDisp$Disparity[1:2])
morpholDisp$sd[3:4] <- sd(morpholDisp$Disparity[3:4])
morpholDisp$sd[5:6] <- sd(morpholDisp$Disparity[5:6])

#Visualize this in ggplot
DispPlot <- ggplot(morpholDisp) + 
  geom_point(aes(y=DispType, x=Disparity, color = Group, group=Group), size = 10, alpha=.8) + 
  scale_color_manual(values = c("skyblue4", "pink3")) + 
  scale_fill_manual(values = c("skyblue4", "pink3")) + 
  labs(x="Morphological Disparity", y=element_blank(), font= "Arial") + 
  geom_errorbar(aes(xmin=Disparity-sd, xmax=Disparity+sd, y=DispType, color = Group, width=0.2)) +
  theme_classic() + 
  theme(legend.position="none", axis.text.y =element_text(colour="black", size = 14), axis.text.x = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = 1) 
DispPlot


############################################
###### DENSITY PLOTS - FIG 1 IN PUB ########
############################################
#Format data -- Note these data are the shape ratios for all traits. No log correction, scaling, etc have been done on this data. 
SR.data.mean_NEW <- SR.data.mean %>%
dplyr::select(-c(SL:Head_W,Pgirdle_to_upr_jaw, X))
row.names(SR.data.mean_NEW) <- SR.data.mean_NEW$Species

SR.data.long <- gather(SR.data.mean_NEW, key="measure", value="value", c("Premax_dentigerous_arm", "AP_of_premax", "Gape_horz", "Jaw_pro", "Buccal_cavity", "LJ_ant.post_position", "LJ_dor.vert_postion", "MN_ant.post_position", "MN_dor.vent_position", "Out_lever", "opening_MA", "closing_MA", "Kinesis", "LPJ_Diag", "LPJ_W", "LPJ_Dwteeth", "AveTeeth", "UPJ_L", "UPJ_W", "UPJ_Dwteeth")) #Format df for facet_wrap

SR.data.long$measure <- gsub("Premax_dentigerous_arm","Premaxillary Dentigerous Arm Length", SR.data.long$measure) 
SR.data.long$measure <- gsub("AP_of_premax","Premaxillary Ascending Process Length", SR.data.long$measure) 
SR.data.long$measure <- gsub("Gape_horz","Gape", SR.data.long$measure) 
SR.data.long$measure <- gsub("Jaw_pro","Oral Jaw Protrusion", SR.data.long$measure) 
SR.data.long$measure <- gsub("Buccal_cavity","Buccal Cavity Length", SR.data.long$measure) 
SR.data.long$measure <- gsub("LJ_ant.post_position","Horizontal Position of Lower Jaw Joint", SR.data.long$measure) 
SR.data.long$measure <- gsub("LJ_dor.vert_postion","Vertical Position of Lower Jaw Joint", SR.data.long$measure) 
SR.data.long$measure <- gsub("MN_ant.post_position","Horizontal Position of Maxilla-Nasal Joint", SR.data.long$measure) 
SR.data.long$measure <- gsub("MN_dor.vent_position","Vertical Position of Maxilla-Nasal Joint", SR.data.long$measure) 
SR.data.long$measure <- gsub("Out_lever","Lower Oral Jaw Length", SR.data.long$measure) 
SR.data.long$measure <- gsub("opening_MA","Opening Mechanical Advantage", SR.data.long$measure) 
SR.data.long$measure <- gsub("closing_MA","Closing Mechanical Advantage", SR.data.long$measure) 
SR.data.long$measure <- gsub("Kinesis","Oral Jaw Kinesis", SR.data.long$measure) 
SR.data.long$measure <- gsub("LPJ_Diag","Lower Pharyngeal Jaw Length", SR.data.long$measure) 
SR.data.long$measure <- gsub("LPJ_W","Lower Pharyngeal Jaw Width", SR.data.long$measure) 
SR.data.long$measure <- gsub("LPJ_Dwteeth","Lower Pharyngeal Jaw Depth", SR.data.long$measure) 
SR.data.long$measure <- gsub("AveTeeth","Lower Pharyngeal Jaw Tooth Diameter", SR.data.long$measure) 
SR.data.long$measure <- gsub("UPJ_L","Upper Pharyngeal Jaw Length", SR.data.long$measure) 
SR.data.long$measure <- gsub("UPJ_W","Upper Pharyngeal Jaw Width", SR.data.long$measure) 
SR.data.long$measure <- gsub("UPJ_Dwteeth","Upper Pharyngeal Jaw Depth", SR.data.long$measure) 

SR.data.long$measure <- as.factor(SR.data.long$measure) #Make all traits a factor 
levels(SR.data.long$measure) #Check levels (i.e., the order of the traits)
SR.data.long$measure <- fct_relevel(SR.data.long$measure, c("Premaxillary Ascending Process Length","Premaxillary Dentigerous Arm Length","Gape","Oral Jaw Protrusion","Oral Jaw Kinesis","Buccal Cavity Length","Lower Oral Jaw Length","Horizontal Position of Lower Jaw Joint","Vertical Position of Lower Jaw Joint","Horizontal Position of Maxilla-Nasal Joint","Vertical Position of Maxilla-Nasal Joint","Opening Mechanical Advantage","Closing Mechanical Advantage","Lower Pharyngeal Jaw Length","Lower Pharyngeal Jaw Width","Lower Pharyngeal Jaw Depth","Lower Pharyngeal Jaw Tooth Diameter","Upper Pharyngeal Jaw Length","Upper Pharyngeal Jaw Width","Upper Pharyngeal Jaw Depth")) #reorder the traits


#Visualize 
trait_density <- ggplot(data=SR.data.long, aes(x=value, colour=Family,  fill=Family)) +
  facet_wrap(~measure, scales="free", labeller = labeller(measure = label_wrap_gen(width = 25)), strip.position = "bottom") +
  #scale_y_continuous(expand=c(0,0)) + #Removes space btwn the data and xaxis
  labs(x="Trait Value", y="Density") +
  geom_density(size=.5, adjust=1.2, alpha = 0.75) + 
  scale_color_manual(values = clr2) + 
  scale_fill_manual(values = clr2) + 
  theme_classic() + 
  theme(legend.position="none", axis.text=element_text(colour="black", size = 8), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .5), axis.ticks.length = unit(.75, "mm"), axis.line=element_line(color = "black", size = .625), aspect.ratio = .65, strip.background=element_blank(),strip.placement="outside", strip.text.x = element_text(size = 10, vjust = 1.25))
trait_density


#Individual plot
Premax_dent_plot <- ggplot(SR.data.mean_NEW, aes(Premax_dentigerous_arm, colour=Family, fill=Family)) +
  geom_density(size=.25, adjust=1.2, alpha = 0.5) + labs(x=paste0("Premaxillary Dentigerous Arm Length", "\n", " "), y = element_blank()) +
  scale_color_manual(values = clr2) + scale_fill_manual(values = clr2) + 
  theme_classic(base_size=20) + 
  theme(legend.position="none", axis.text =element_text(colour="black", size = 6.5), axis.title=element_text(colour="black", size = 10), axis.ticks=element_line(colour="black", size= .25), axis.ticks.length = unit(1, "mm"), axis.line=element_line(color = "black", size = .45), aspect.ratio = .75)
```




```{r RATES OF EVOLUTION, warning=FALSE}   

# <--------------- FORMAT DATA ---------------> #
disp_Combo3Fin <- prepca %>% 
  mutate_all(funs(as.numeric(.))) 
disp_Combo3Fin.2 <- as.matrix(disp_Combo3Fin)

disp_Combo3Fin_OJ <- ComboFin2_OJ %>% 
  mutate_all(funs(as.numeric(.))) 
disp_Combo3Fin_OJ.2 <- as.matrix(disp_Combo3Fin_OJ)

disp_Combo3Fin_PJ <- ComboFin2_PJ %>% 
  mutate_all(funs(as.numeric(.))) 
disp_Combo3Fin_PJ.2 <- as.matrix(disp_Combo3Fin_PJ)

GroupFin <- as.factor(ComboFin2_new$Fam2) #Cichs are group PHGY and cents are group non-PHGY
names(GroupFin) <- ComboFin2_new$Species3


# <--------------- RUN MULTIVARIATE RATE ANALYSES in GEOMORPH---------------> #
#Results reported in Table 1
(OverallRate <- compare.evol.rates(disp_Combo3Fin.2, ComboTreeFin2, gp=GroupFin, iter = 10000, method = "simulation", print.progress = TRUE))

(OJRates <- compare.evol.rates(disp_Combo3Fin_OJ.2, ComboTreeFin2, gp=GroupFin, iter = 10000, method = "simulation", print.progress = TRUE))

(PJRates <- compare.evol.rates(disp_Combo3Fin_PJ.2, ComboTreeFin2, gp=GroupFin, iter = 10000, method = "simulation", print.progress = TRUE))


# <--------------- RUN UNIVARIATE RATE ANALYSES in GEOMORPH---------------> #
#Results reported in Fig 1
#Set up data frame to catch results using "groups" and "traitvec" variables that we made in disparity section
AllRates_inter <- as.data.frame(matrix(ncol = ncol(disp_Combo3Fin.2) , nrow = length(groups)+3))
colnames(AllRates_inter) <- c(trait_vec)
rownames(AllRates_inter) <- c(groups,"rate ratio","pairwise pvalue","effect size")

#sequential_disparity
print("starting to run rate")

for(i in 1:ncol(disp_Combo3Fin.2)){
  
  print(paste("Now running rate for", trait_vec[i]))
  
  run_rate_inter <- compare.evol.rates(disp_Combo3Fin.2[,i], ComboTreeFin2, gp=GroupFin, iter = 10000, method = "simulation", print.progress = F)
  
  AllRates_inter[1,i] <- run_rate_inter$sigma.d.gp[1]
  AllRates_inter[2,i] <- run_rate_inter$sigma.d.gp[2]
  AllRates_inter[3,i] <- run_rate_inter$sigma.d.ratio
  AllRates_inter[4,i] <- run_rate_inter$pairwise.pvalue
  AllRates_inter[5,i] <- run_rate_inter$Z
}


# <--------------- VISUALIZE RATE RESULTS ---------------> #
#Not used in pub
#Retrieve variance info + format in df
morpholRates <- as.data.frame(c(OverallRate$sigma.d.gp,OJRates$sigma.d.gp,PJRates$sigma.d.gp))
names(morpholRates) [1]  <- "Rate"
morpholRates$Group <- as.factor(c("Cichlidae", "Centrarchidae", "Cichlidae", "Centrarchidae", "Cichlidae", "Centrarchidae"))
morpholRates$RateType<- as.factor(c("Overall", "Overall", "Oral Jaws", "Oral Jaws", "Pharyngeal Jaws", "Pharyngeal Jaws"))
morpholRates$RateType <- fct_relevel(morpholRates$RateType, c("Pharyngeal Jaws", "Oral Jaws", "Overall"))

#Compute std dev and add it to the df
morpholRates$sd[1:2] <- sd(morpholRates$Rate[1:2])
morpholRates$sd[3:4] <- sd(morpholRates$Rate[3:4])
morpholRates$sd[5:6] <- sd(morpholRates$Rate[5:6])

#Visualize these results 
ratePlot <- ggplot(morpholRates) + 
  geom_point(aes(y=RateType, x=Rate, color = Group, group=Group), size = 10, alpha=.8) + 
  scale_color_manual(values = c("skyblue4", "pink3")) + 
  scale_fill_manual(values = c("skyblue4", "pink3")) + 
  labs(x="Rates of Morphological Evolution", y=element_blank(), font= "Arial") + 
  geom_errorbar(aes(xmin=Rate-sd, xmax=Rate+sd, y=RateType, color = Group, width=0.2)) +
  theme_classic() + 
  theme(legend.position="none", axis.text.y =element_text(colour="black", size = 14), axis.text.x = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = 1) 
 
```




```{r INTEGRATION + MODULARITY, warning=FALSE}   

############################################
###### BETWEEN JAW UNIT INTEGRATION ########
############################################
# <--------------- FORMAT DATA ---------------> #
#Create just centrarchids and cichlid trees 
centRows3 <- Species3[86:115]
cent_tree3 <- drop.tip(ComboTreeFin2, setdiff(ComboTreeFin2$tip.label,centRows3))
cent_tree3$tip.label
cent_treeFin3 <- ladderize(cent_tree3, right=FALSE)

cichRows3 <- Species3[1:85]
cich_tree3 <- drop.tip(ComboTreeFin2, setdiff(ComboTreeFin2$tip.label,cichRows3))
cich_tree3$tip.label
cich_treeFin3 <- ladderize(cich_tree3, right=FALSE)


#Create matrices for each jaw module in cichlids and centrarchids (4 df total). Standardize/scale these values so they are z-scores -- this makes sure that the values are in commensurate units and scales. Standardize the groups separately as they may be on different scales, areas of morphospace, etc. 
integ_mod_Cich <- prepca[c(1:85),] %>% 
  mutate_all(funs(scale(.)))

integ_mod_Cent <- prepca[c(86:115), ] %>% 
  mutate_all(funs(scale(.)))

integ_mod_mat <- rbind(integ_mod_Cich,integ_mod_Cent)
integ_mod_mat2 <- as.matrix(integ_mod_mat)

Inter_OJCich6 <- as.matrix(integ_mod_mat[c(1:85),c(1:13)])
Inter_PJCich6 <- as.matrix(integ_mod_mat[c(1:85),c(14:20)])

Inter_OJCent6 <- as.matrix(integ_mod_mat[c(86:115),c(1:13)])
Inter_PJCent6 <- as.matrix(integ_mod_mat[c(86:115),c(14:20)])


#Check the structure and class of the modules. Structure should be characters for all the row and column names and the numbers should be numeric. Class should be matrix or array
str(Inter_PJCich6)
class(Inter_PJCich6)


# <--------------- RUN BETWEEN-JAW-UNIT INTEGRATION ANALYSES ---------------> #
#Results reported in Table 2
#Calculate integration for centrarchids and cichlids -- Function quantifies the degree of phylogenetic morphological covariation between two or more sets of Procrustes shape variables using partial least squares. This tells us how independent the OJ and PJ biomechanical units are from each other in each family. 
CentInter3 <- phylo.integration(Inter_OJCent6, Inter_PJCent6 ,partition.gp = NULL, phy=cent_treeFin3, iter = 9999, seed = NULL, print.progress = TRUE)
CentInter3

CichInter3 <- phylo.integration(Inter_OJCich6, Inter_PJCich6 ,partition.gp = NULL, phy=cich_treeFin3, iter = 9999, seed = NULL, print.progress = TRUE)
CichInter3

#Compare the amount of independence found among the jaw modules in cichlids vs centrarchids. This function determines whether the two r-PLS values are significantly different. PLS — "This method rotates each multivariate dataset to the axis of most covariance with the other set of landmark coordinates and tests for a correlation between these two morphological axes.  (From ronco and salzburger) So this is very very similar to a PC axis but its the PC1 axis of one block of traits vs that of the other block of traits" -- From Ronco & Salzburger, 2020
compare.pls(CentInter3,CichInter3)


# <--------------- VISUALIZE THE RESULTS ---------------> #
#Format data
PLS_CichScores <- as.data.frame(CichInter3$XScores[ ,1])
names(PLS_CichScores)[1] <- "X"
PLS_CichScores$Y <- CichInter3$YScores[ ,1]

PLS_CentScores <- as.data.frame(CentInter3$XScores[ ,1])
names(PLS_CentScores)[1] <- "X"
PLS_CentScores$Y <- CentInter3$YScores[ ,1]

Patterns_PLS <- rbind(PLS_CichScores, PLS_CentScores)
Patterns_PLS$Fam <- ComboFin2_new$Fam2


#Cichlid results
CichPattents_PLS <- ggplot(data = PLS_CichScores, aes(x = X, y = Y, col = "rosybrown")) + 
  geom_point( size = 5) +
  scale_colour_manual(values = "rosybrown") +
  xlab("Block 1 ~ PLS 1") +
  ylab("Block 2 ~ PLS 2") +
  theme_classic() +
  theme(legend.position="none") + 
  geom_smooth(method = lm , se = FALSE, col = "black", size=2) #+
#geom_text(label=(ComboFin2$Species3[1:85]), size=1.5, hjust=0, vjust=0)
CichPattents_PLS


#Centrarchid results
CentPattents_PLS <- ggplot(data = PLS_CentScores, aes(x = X, y = Y, col = "skyblue4")) + 
  geom_point( size = 5) +
  scale_colour_manual(values = "skyblue4") +
  xlab("Block 1 ~ PLS 1") +
  ylab("Block 2 ~ PLS 2") +
  theme_classic() +
  theme(legend.position="none") + 
  geom_smooth(method = lm, se = FALSE, col = "black", size=2) #+
#geom_text(label=(ComboFin2$Species3[86:115]), size=1.5, hjust=0, vjust=0)
CentPattents_PLS


#Combined reuslts -- Figure 5
Pattents_PLSplot <- ggplot(data = Patterns_PLS, aes(x = X, y = Y, col = Fam2)) + 
  geom_point(size = 3, alpha=.8) + 
  geom_smooth(method = lm, se = FALSE, size=1.5) +
  labs(x= "PLS Block 1 ~ Oral Jaw Scores", y= "PLS Block 2 ~ Pharyngeal Jaw Scores", font="Arial") + 
  scale_colour_manual(values = clr2) + 
  scale_fill_manual(values = clr2) + 
  theme_classic() + 
  theme(legend.position="none", axis.text = element_text(colour="black", size = 10), axis.title=element_text(colour="black", size = 14), axis.ticks=element_line(colour="black", size= .75), axis.ticks.length = unit(1.5, "mm"), axis.line=element_line(color = "black", size = .75), aspect.ratio = 1)
  #+
#geom_text(label=(ComboFin2$Species3[86:115]), size=1.5, hjust=0, vjust=0)



##########################
###### MODULARITY ########
##########################
#Are the oral jaws and pharyngeal jaws true modules in cichlids and centrarchids? In other words, is the integration within each jaw module stronger than the integration between the modules??? Modularity describes patterns where trait correlations are high within subsets of variables, termed modules, and where trait correlations across modules are comparatively weaker (Adams & Collyer, 2019).

#The function below quantifies the degree of modularity between two or more hypothesized modules of Procrustes shape variables in a phylogenetic context and compares this to patterns found by randomly assigning landmarks into subsets. The degree of modularity is quantified using the CR coefficient (Adams 2016; calculated as C(S1S2)/[C(S1)*C(S2)]). The CR coefficient for the observed modular hypothesis is then compared to a distribution of values obtained by randomly assigning landmarks into subsets, with the restriction that the number of landmarks in each subset is identical to that observed in each of the original partitions. A significant modular signal is found when the observed CR coefficient is small relative to this distribution (see Adams 2016). Such a result implies that there is significantly greater independence among modules than is expected under the null hypothesis of random associations of variables (neither modular nor integrated structure). This result is consistent with the identification of significant modular structure in the data. In addition, a multivariate effect size describing the strength of the effect is estimated from the empirically-generated sampling distribution (see details in Adams and Collyer 2019). A histogram of coefficients obtained via resampling is presented, with the observed value designated by an arrow in the plot.

#Empirically, the CR coefficient has a lower limit of 0, an expected value of 1 when there is no modularity, and values greater than 1 are also possible if covariances between module variables exceed covariances within modules. Thus, CR values that are closer to 0 describe datasets with a relatively greater degree of modularity, while CR values closer to 1 describe datasets with less modular signal. To evaluate CR statistically, a resampling procedure that randomly assigns variables to modules is used, and in each iteration a CRrand is obtained under the null hypothesis of no modularity. The observed value, CRobs, is then compared to the empirically-generated sampling distribution of CRrand values to evaluate significance (for details see: Adams 2016). This will give us a Zcr or effect size score. One should recognize that more negative values of ZCR represent greater modular signal, because when modular signal is present, CR will be less than μˆ .

# <--------------- FORMAT DATA ---------------> #
jaw_mod3 <- rep('OJ',20); jaw_mod3[14:20]<-'PJ'

Mod_Cent4 <- as.matrix(integration_ANOVA_mat2[c(86:115), ])
Mod_Cich4 <- as.matrix(integration_ANOVA_mat2[c(1:85), ])


# <--------------- RUN MODULARITY ANALYSES ---------------> #
#Results reported in Table 2
centMod3 <- phylo.modularity(Mod_Cent4,jaw_mod3,phy=cent_treeFin3,iter = 10000, CI = FALSE, seed = NULL,print.progress = TRUE)
centMod3
plot(centMod3)

cichMod3 <- phylo.modularity(Mod_Cich4,jaw_mod3, phy=cich_treeFin3,iter = 10000, CI = FALSE, seed = NULL,print.progress = TRUE)
cichMod3
plot(cichMod3)


#Determine whether the strength of modularity is significantly different across datasets. Here we will compare the strength of modular signal (CR) across cichlids and centrarchids. We will use a two-sample test, which calculates the effect size of the diverence between pairs of modular effect sizes.  
modStren3 <- compare.CR(centMod3, cichMod3, CR.null = FALSE, two.tailed = TRUE)
modStren3$pairwise.z #A matrix of pairwise, two-sample z scores between all pairs of effect sizes
modStren3$pairwise.P #A matrix of corresponding P-values



##########################################
###### WITIN JAW UNIT INTEGRATION ########
##########################################
####Integration within each jaw unit. Basically, one can estimate the magnitude (strength) of integration in a set of traits using the relative eigenvalue index (Vrel). Then, we convert this to an effect size (Z-score), and then use these to compare the strength of integration between samples. ---> "The function quantifies the strength of morphological integration in a set of variables. Here the set of traits are treated as a single unit, and the overall degree of covariation in them is quantified using the relative eigenvalue index: Vrel (Pavlicev et al. 2009). Following Conaway and Adams (2022), only the non-trivial dimensions of variation are used in the calculation of Vrel. The measure is then converted to an effect size (Z-score), based on the procedures in Conaway and Adams (2022). These may be used in subsequent comparisons of the strength of integration across datasets. Input for the analysis may be a 3D array of Procrustes coordinates, of a matrix of variables. If the observations are species related by a phylogeny, the phylogeny may also be included." From the help page!

# <--------------- RUN INTEGRATION ANALYSES ---------------> #
#Results reported in Table 2
#Oral jaw integration
(CentInter_OJ <- integration.Vrel(Inter_OJCent6, phy=cent_treeFin3))
(CichInter_OJ <- integration.Vrel(Inter_OJCich6, phy=cich_treeFin3))

#Pharyngeal jaw integration
(CentInter_PJ <- integration.Vrel(Inter_PJCent6, phy=cent_treeFin3))
(CichInter_PJ <- integration.Vrel(Inter_PJCich6, phy=cich_treeFin3))

#Compare results
OJ_comparison <- compare.ZVrel(CentInter_OJ,CichInter_OJ)
OJ_comparison

PJ_comparison <- compare.ZVrel(CentInter_PJ,CichInter_PJ)
summary(PJ_comparison)
```




```{r COEVOLUTIONARY PATTERNS, warning=FALSE}   
##### Visualize correlation matrices, pairwise correlations between trait values and differences in Cichlids and Centrarchids####
#Code below is modified from Dean Adams (emails and gitHub repository) --> A is a 2D matrix, rows = species (with rownames as species names), and phy is the phylogeny. NOTE: MAKE SURE YOUR TREE ONLY INCLUDES SPECIES IN THE DATASET

# <--------------- CREATE RELEVANT FUNCTIONS ---------------> #
#Underlying functions from github repository
phy.sim.mat <- function(phy) {
  N <- length(phy$tip.label)
  n <- nrow(phy$edge)
  m <- matrix(0, N, n)
  edg <- phy$edge.length
  idx <- phy$edge[, 2]
  anc <- phy$edge[, 1]
  tips <- which(idx <= N)
  non.tips <- which(idx > N)
  for(i in 1:length(tips)) m[idx[tips[i]], tips[i]] <- sqrt(edg[tips[i]])
  for(i in 1:length(non.tips)) {
    x <- idx[non.tips[i]]
    anc.i <- which(anc == x)
    edg.i <- idx[anc.i]
    if(any(edg.i > N)) {
      edg.i <- as.list(edg.i)
      while(any(edg.i > N)) {
        edg.i <- lapply(1:length(edg.i), function(j){
          edg.i.j <- edg.i[[j]]
          if(edg.i.j > N) {
            idx[which(anc == edg.i.j)]
          } else edg.i.j <- edg.i.j
        })
        edg.i <- unlist(edg.i)
      }
    }
    m[edg.i, which(idx == x)] <- sqrt(edg[which(idx == x)])
  }
  rownames(m) <- phy$tip.label
  m
}

fast.phy.vcv <- function (phy) tcrossprod(phy.sim.mat(phy)) ## same as vcv.phylo but without options, in order to not use ape

fast.solve <- function(x) { 
  res <- try(solve(x), silent = TRUE)
  if(inherits(res, "try-error")) res <- fast.ginv(x)
  return(res)
} # same as solve, but without traps (faster)

fast.ginv <- function(X, tol = sqrt(.Machine$double.eps)){
  k <- ncol(X)
  if(inherits(X, "matrix")) {
    Xsvd <- La.svd(X, k, k)
    Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
    rtu <-((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE]))
    v <-t(Xsvd$vt)[, Positive, drop = FALSE]
  } else {
    Xsvd <- svd(X, k, k)
    Positive <- Xsvd$d > max(tol * Xsvd$d[1L], 0)
    rtu <-((1/Xsvd$d[Positive]) * t(Xsvd$u[, Positive, drop = FALSE]))
    v <-Xsvd$v[, Positive, drop = FALSE]
  }
  
  v %*% rtu
} # same as ginv, but without traps (faster)

phylo.mat<-function(x,phy){
  C<-fast.phy.vcv(phy)
  C<-C[rownames(x),rownames(x)]
  invC <-fast.solve(C)
  eigC <- eigen(C)
  if(any(Im(eigC$values)==0)) eigC$values<-Re(eigC$values)
  if(any(Im(eigC$vectors)==0)) eigC$vectors<-Re(eigC$vectors)
  lambda <- zapsmall(abs(Re(eigC$values)))
  if(any(lambda == 0)){
    warning("Singular phylogenetic covariance matrix. Proceed with caution")
  }
  D.mat <- fast.solve(eigC$vectors%*% diag(sqrt(abs(eigC$values))) %*% t(eigC$vectors))
  rownames(D.mat) <- colnames(D.mat) <- colnames(C)
  list(invC = invC, D.mat = D.mat,C = C)
}


# <--------------- FORMAT DATA ---------------> #
#Cent
disp_Combo_Cent <- as.matrix(disp_Combo3Fin.2[86:115, ])
n <- nrow(disp_Combo_Cent)
phy.parts <- phylo.mat(disp_Combo_Cent, cent_treeFin3)
Ptrans <- phy.parts$D.mat %*% (diag(n) - matrix(1, n) %*%
                                 crossprod(matrix(1, n), phy.parts$invC)/sum(phy.parts$invC))
Matx_Cent <- Ptrans %*% disp_Combo_Cent

#Cich
disp_Combo_Cich <- as.matrix(disp_Combo3Fin.2[1:85, ])
n2 <- nrow(disp_Combo_Cich)
phy.parts2 <- phylo.mat(disp_Combo_Cich, cich_treeFin3)
Ptrans2 <- phy.parts2$D.mat %*% (diag(n2) - matrix(1, n2) %*%
                                 crossprod(matrix(1, n2), phy.parts2$invC)/sum(phy.parts2$invC))
Matx_Cich <- Ptrans2 %*% disp_Combo_Cich


# <--------------- CALCULATE + VISUALIZE CORRELATION + COVARIATION MATRICES ---------------> #
#Calculate the correlation + covariance matrices
CorrMatx_Cents <- cor(Matx_Cent)
CorrMatx_Cichs <- cor(Matx_Cich)
CovMatx_Cents <- cov(Matx_Cent)
CovMatx_Cichs <- cov(Matx_Cich)


#Visualize correlation matrices -- Fig 4a
CorrCols <- colorRampPalette(c("dark blue","light blue","gray", "white", "orange", "red", "dark red"))(7)

Cich_Corr_plot <- corrplot(CorrMatx_Cichs[], method = "circle", type = "upper", col = CorrCols, number.cex = .5, tl.col = "black", cl.cex = .8, cl.ratio = .1, cl.pos = "r")
Cich_Corr_plot 

Cent_Corr_plot <- corrplot(CorrMatx_Cents[], method = "circle", type = "lower", col = CorrCols, number.cex = .5, tl.col = "black", cl.cex = .8, cl.ratio = .1, cl.pos = "r")
Cent_Corr_plot 


#Visualize the abs. value difference between the two matrices -- Fig 4b
CorrCols2 <- colorRampPalette(c("skyblue4","skyblue1","white","rosybrown1","rosybrown4"))(7)

ABSCorrDiff <- as.matrix(abs(CorrMatx_Cichs) - abs(CorrMatx_Cents))
CorrDiff2 <- as.data.frame(ABSCorrDiff)
Diff_Corr_plot <- corrplot(ABSCorrDiff[], method = "circle", type = "lower", col = CorrCols2, is.corr = FALSE, tl.col = "black", cl.cex = .8, cl.ratio = .1, cl.pos = "r")
```