#################################
### SCRIPT BENDER ###############
################################

###########################################
# FUNCTIONAL DIVERSITY OF CHAETODONTIDAE #
##########################################


library(FD)
library(vegan)
library(cluster)
library(reshape2)
library (plyr)

#Importing datasets ########

#set you working directory below

setwd("/Users/marianabender/Desktop/Func_chaetodon")

data <- read.csv("chaet_dist.csv", header=T, sep=";")
traits <- read.table("trait.txt",T)

########################
#Take a look at data ###
########################
head (data)
head(data[1:3, 1:3])
head(traits)

########################
#merge two datasets ###
########################

data_chaet = merge(data, traits, by="Genus_and_species")


###### Calculating functional indices for Chaetodontidae  ###################
############# Atlantic vs. Indo-Pacific species  ############################
#############################################################################

# first we will apply a dissimilarity analysis using the gower distance
# applyied to functional traits

data_chaet=na.omit(data_chaet) #loose species for which there are NAs

trait_chaet = data_chaet[,c("Size","Diet_2012",
                            "Home_range","Activity",
                            "Schooling","Level")]


#now, run a gower distance for all Chaetodontidae species

gower_all = daisy (trait_chaet, metric="gower")
pcoa_all = cmdscale(gower_all, k=4)

#paste pcoa axes in data set

data_chaet = cbind(data_chaet, pcoa_all)


############# Plot functional volumes for Atlantic vs. Indo-Pacific species  ############################

##
par(mfrow=c(3,2))
plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Indo-Pacific") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", alpha=70, col="gold", bg="black", show.groups=TRUE)
ordihull(pcoa_all, groups=data_chaet$Indian>=1 | data_chaet$Pacific>=1, lty="dashed", draw="polygon", alpha=100, col="magenta", show.groups=TRUE)


plot (pcoa_all[,3],pcoa_all[,4], xlab="PC3", ylab="PC4", cex=0.9, xlim=c(-0.7,0.5), ylim=c(-0.4,0.7),pch=21, bg="white") 
ordihull(pcoa_all[,3:4], groups=data_chaet$Size_class>=1, lty="dashed", draw="polygon", alpha=70, col="gold", bg="black", show.groups=TRUE)
ordihull(pcoa_all[,3:4], groups=data_chaet$Indian>=1 | data_chaet$Pacific>=1, lty="dashed", draw="polygon", alpha=100, col="magenta", show.groups=TRUE)

##
plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Atlantic") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", alpha=70, col="gold", bg="black", show.groups=TRUE)
ordihull(pcoa_all, groups=data_chaet$Atlantic>=1, lty="dashed", draw="polygon", alpha=200, col="darkblue", show.groups=TRUE)

plot (pcoa_all[,3],pcoa_all[,4], xlab="PC3", ylab="PC4", cex=0.9, xlim=c(-0.7,0.5), ylim=c(-0.4,0.7),pch=21, bg="white") 
ordihull(pcoa_all[,3:4], groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", alpha=70, col="gold", bg="black", show.groups=TRUE)
ordihull(pcoa_all[,3:4], groups=data_chaet$Atlantic>=1, lty="dashed", draw="polygon", alpha=200, col="darkblue", show.groups=TRUE)

data_chaet[data_chaet$`1`>=0.2031 & data_chaet$`1`<=0.4 &
             data_chaet$`2`<=-0.1 & data_chaet$`2`>=-0.3, ]
## Now, check the volume filled by Indo-Pacific vs. Atlantic Butterflyfishes

# We can also explore how Butterflyfishes from the Tropical Eastern Pacific occupy the functional space

#par(mfrow=c(1,2))
plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="TEP") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", alpha=70, col="gold", bg="black", show.groups=TRUE)
ordihull(pcoa_all, groups=data_chaet$TEP>=1, lty="dashed", draw="polygon", alpha=200, col="palegreen4", show.groups=TRUE)

plot (pcoa_all[,3],pcoa_all[,4], xlab="PC3", ylab="PC4", cex=0.9, xlim=c(-0.7,0.5), ylim=c(-0.4,0.7),pch=21, bg="white") 
ordihull(pcoa_all[,3:4], groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", alpha=70, col="gold", bg="black", show.groups=TRUE)
ordihull(pcoa_all[,3:4], groups=data_chaet$TEP>=1, lty="dashed", draw="polygon", alpha=200, col="palegreen4", show.groups=TRUE)

## Using Ordiellipse function we can explore how traits appear in the functional space

par(mfrow=c(2,3))
plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Sessile invertivores") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", bg="black", col="lightgray", show.groups=TRUE)
ordiellipse(pcoa_all, groups=data_chaet$Diet_2012=="IS", lty=1, col="pink", draw="polygon", show.groups=TRUE)

plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Mobile invertivores") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", bg="black", col="lightgray",show.groups=TRUE)
ordiellipse(pcoa_all, groups=data_chaet$Diet_2012=="IM", lty=1, col="darkgreen", draw="polygon", show.groups=TRUE)

plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Planktivores") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", bg="black", col="lightgray",show.groups=TRUE)
ordiellipse(pcoa_all, groups=data_chaet$Diet_2012=="PK", lty=1, col="magenta", draw="polygon", show.groups=TRUE)

plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Solitary") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", bg="black", col="lightgray",show.groups=TRUE)
ordiellipse(pcoa_all, groups=data_chaet$Schooling=="S", lty=1, col="orange", draw="polygon", show.groups=TRUE)

plot (pcoa_all[,1],pcoa_all[,2], xlab="PC1", ylab="PC2", cex=0.9, xlim=c(-0.4,0.7), ylim=c(-0.4,0.7),pch=21, bg="white", main="Pairing") 
ordihull(pcoa_all, groups=data_chaet$Size_class>=1, lty="solid", draw="polygon", bg="black", col="lightgray",show.groups=TRUE)
ordiellipse(pcoa_all, groups=data_chaet$Schooling=="P", lty=1, col="blue", draw="polygon", show.groups=TRUE)

### check which set of traits are potentially missing from the Atlantic and TEP


############################################
# Calculate Functional Indices  ############
############################################

#We will apply dbFD index to compute functional diversity indices for fish communities


##########################################
# Analysis at the realm scale ############
##########################################

#trait data
trait_chaet = data_chaet[,c("Size","Diet_2012","Home_range","Activity","Schooling","Level")]
trait_chaet$Size <- as.numeric(as.character(trait_chaet$Size))

#adjust data row names
row.names(trait_chaet)=paste(data_chaet$Genus_and_species)

#apply gower distance
dist = gowdis(trait_chaet, ord="podani") #ordinal variables will be ranks

#community matrix 1
comm = data_chaet[170:173]
row.names(comm)=paste(data_chaet$Genus_and_species)

## Run dbFD
chaet_func = dbFD(dist, t(comm), w.abun = FALSE,
                  ord = c("podani", "metric"),
                  corr = "cailliez",
                  calc.FRic = TRUE, stand.FRic = TRUE)

##########################################
# Analysis at the site scale #############
##########################################

#community matrix 2

comm2 = data_chaet[1:169]
rownames(comm2) <- comm2$Genus_and_species
comm2 <- comm2[,-1]
comm3 = comm2[,colSums(comm2) >= 5] #select sites for which there are >5 Chaetodontidae spp
comm4 = comm3[rowSums(comm3) >= 1,] #eliminate sites with 0 spp

#trait data
trait_chaet = data_chaet[,c("Genus_and_species","Size","Diet_2012","Home_range","Activity","Schooling","Level")]
trait_chaet$Size <- as.numeric(as.character(trait_chaet$Size))

comm4$Genus_and_species <- rownames (comm4)
rownames (comm4) <- NULL
comm4 = merge(comm4, trait_chaet, by= "Genus_and_species")

#extract community data
comm_site = comm4[,2:126]
row.names(comm_site) = paste(comm4$Genus_and_species)

#extract traits and run gower
comm_traits = comm4[,127:132]
row.names(comm_traits) = paste(comm4$Genus_and_species)

#apply gower distance
dist2 = gowdis(comm_traits, ord="podani")

#Genus_and_Species(comm3)=paste(data_chaet$Genus_and_species)

## run dbFD
chaet_func2 = dbFD(dist2, t(comm_site), w.abun = FALSE, ord = c("podani", "metric"),
                   corr = "cailliez", calc.FRic = TRUE, stand.FRic = TRUE)

## check out chaet_func2 object
chaet_func2

##################################################
# Plot functional indices for Chaets #############
##################################################

# Now we will explore how the Functional Indices vary according to local species 
# richness

par(mfrow=c(2,2))

#Functional richness
plot(chaet_func2$nbsp, chaet_func2$FRic, ylab="Functional Richness", xlab="Species Richness",
      col="black", bg="blue", pch=21, cex=1.2)

#Now lets take a look at functional evenness
plot(chaet_func2$nbsp, chaet_func2$FEve, ylab="Functional Evenness", ylim=c(0.5,1.0), xlab="Species Richness",
     col="black", bg="orange", pch=21, cex=1.2)

#Now lets take a look at functional divergence
plot(chaet_func2$nbsp, chaet_func2$FDiv, ylab="Functional Divergence", ylim=c(0.5,1.0), xlab="Species Richness",
     col="black", bg="red", pch=21, cex=1.2)


##################################################