####
# Figure 2 - Benthic cover at different depths and protection
####

# set up
rm(list=ls())

library(reshape)
library(vegan)

# Reading in data

df1<- read.csv("ESM 2 Benthic data.csv",
               header=TRUE, colClasses='factor')


##
# Tallying up how many points of each benthic type are recored on each transect
##

# Pasting the site and transect together to create a unqiue
# identifier for each transect

df1$Location <- paste (df1$Site, df1$Depth, df1$Transect, sep="_")

# Generating frequencies of each location
df2 <- table(c(df1$Location), df1$Category)
df2 <- as.data.frame(df2)
colnames(df2) <- c("Location", "Category", "Frequency")


##
# Calculate percentage cover for each catagory
##

df2$Percentage <- NA

for (i in unique(df2$Location))
{
  
  # Calculating total data points for this location
  temp <- df2[which(df2$Location == i), ]
  total.points <- sum(temp$Frequency)
  
  for (j in unique(df2$Category))
  {
    
    # Calculating the percentage cover of each category
    row.no <- which(df2$Location == i & df2$Category == j)
    df2$Percentage[row.no] <- (df2[row.no, "Frequency"] / total.points) * 100
    
  } # closes the Category loop
} # closes the Location loop


# seperating the location column into site, depth and transect
df2$Site <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 1)
df2$Depth <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 2)
df2$Transect <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 3)


# removing the now redundant location column (first column)
df2 <- subset(df2, select= -1)


# Function to calculate benthic cover of any catagory

benthic_cover <- function(code, input.data)
{
  benthic <- input.data[which(input.data$Category==code), ]
  
  # Results dataframe to hold mean cover for each site and depth and then SE
  
  benthic_summary <- data.frame(sort(unique(as.numeric(benthic$Depth))))
  colnames(benthic_summary) <- "Depth"
  benthic_summary$Site <- rep(NA, nrow(benthic_summary))
  benthic_summary$Category <- rep(NA, nrow(benthic_summary))
  benthic_summary$Mean <- rep(NA, nrow(benthic_summary))
  benthic_summary$SE <- rep(NA, nrow(benthic_summary))
  
  # Collapsing the dataframe back to having no rows so
  # rows can be populated as results are generated
  benthic_summary <- benthic_summary[0, ]
  
  # Calculating mean cover at each site and depth and then SE
  
  # looping through sites
  for (site in unique(as.character(benthic$Site)))
  {
    temp_site <- benthic[which(benthic$Site == site), ]
    
    # looping through depths
    for (depth in sort(unique(as.numeric(benthic$Depth))))
    {
      temp <- temp_site[which(temp_site$Depth==depth), ]
      
      # results row number in the output dataframe
      row <- nrow(benthic_summary) + 1
      
      benthic_summary[row, c("Site", "Depth", "Category")] <- c(site, depth, code)
      
      #mean
      benthic_summary[row, 'Mean'] <- sum(temp$Percentage)/length(temp$Percentage)
      
      #se
      benthic_summary[row, "SE"] <- sqrt(var(temp$Percentage)/length(temp$Percentage))
      
    } # closes the depth in loop
  } # closes the site loop
  
  return (benthic_summary)
} # closes benthic_cover function



##
# Calculating % cover of different types
##

# Hard Coral
df3 <- benthic_cover(code="HC", input.data=df2)

# Sponge
df3 <- rbind(df3, benthic_cover(code="SPON", input.data=df2))

# Gorgonians
df3 <- rbind(df3, benthic_cover(code="GORG", input.data=df2))

# Macroalgae Fleshy
df3 <- rbind(df3, benthic_cover(code="MAF", input.data=df2))

# Macroalgae Calcarious
df3 <- rbind(df3, benthic_cover(code="MAC", input.data=df2))

# Black Coral
df3 <- rbind(df3, benthic_cover(code="BC", input.data=df2))

# Sand, Rubble and Pavement
df3 <- rbind(df3, benthic_cover(code="OSUB", input.data=df2))


##
# Adding Protection Status to Depth Column to seperate out results based on 
##

protected <- c("CO", "HE", "PJ", "PT", "SR")
unprotected <- c("PU", "TT", "VB")

for (i in 1:nrow(df3))
{
  
  # Protected sites
  if (df3$Site[i] %in% protected)
  {
    df3$Depth[i] <- paste(df3$Depth[i], "PRO", sep="_")
  }
  
  # Unprotected sites
  if (df3$Site[i] %in% unprotected)
  {
    df3$Depth[i] <- paste(df3$Depth[i], "UNP", sep="_")
  }
  
} # closes the i in 1:nrow(df3) loop



##
# Calculating summary statistics of each depth and code
##

# results table
col <- which(colnames(df3)=="Site")
df4 <- subset(df3, select= -col)
df4 <- df4[0,]

for (code in unique(df3$Category))
{
  for (depth in unique(df3$Depth))
  {
    temp <- df3[which(df3$Depth == depth & df3$Category == code), ]
    
    # results row number in the output dataframe
    row <- nrow(df4) + 1
    
    df4[row, c("Depth", "Category")] <- c(depth, code)
    
    #mean
    df4[row, 'Mean'] <- sum(temp$Mean)/length(temp$Mean)
    
    #se
    df4[row, "SE"] <- sqrt(var(temp$Mean)/length(temp$Mean))
    
  } # closes the depth loop
} # closes the benthic code loop


##
# Seperating the Depth and Protection Status
##

# seperating the location column into site, depth and transect for df
df3$Protection <- sapply(strsplit(as.character(df3$Depth), '_'), "[[", 2)
df3$Depth <- sapply(strsplit(as.character(df3$Depth), '_'), "[[", 1)


# seperating the location column into site, depth and transect for df4
df4$Protection <- sapply(strsplit(as.character(df4$Depth), '_'), "[[", 2)
df4$Depth <- sapply(strsplit(as.character(df4$Depth), '_'), "[[", 1)


##
# Stats Function - Wilcox test between benthic groups
##

difference <- function(code, depth, df)
{
  cover1 <- df[which(df$Category==code & df$Depth == depth & df$Protection=='PRO'), ]
  cover2 <- df[which(df$Category==code & df$Depth == depth & df$Protection=='UNP'), ]
  
  return(wilcox.test(cover1$Mean, cover2$Mean))
  
} # closes the difference function

# Example of stats for 15m Sponges
# difference("SPON", 15, df3)

####
# Plotting
####

# splitting the dataframe into Shallow and Deep for plotting

shal <- df4[which(df4$Depth=="15"), ]
deep <- df4[which(df4$Depth=="55"), ]


###
# Plot set up
##

tiff("Fig 2  - Benthic.tiff", width = 6, height = 6, units = 'in', res = 300)


par(oma=c(5, 2, 0, 0))
par(mfrow=c(2, 1))
par(mar=c(2, 2, 1, 0))
par(las=2)


##
# Plotting shallows
##

barplot(matrix(shal$Mean, nc=length(unique(shal$Category))), beside=TRUE,
        axes=FALSE, ylim=c(0, 80), col=c('grey100', 'grey66'))

axis(1, at=c(2, 5, 8, 11, 14, 17, 20), labels=rep("", 7), cex=0.3, lwd.ticks=FALSE)

axis(2)

# Adding error bars

location <- c(1.5, 2.5, 4.5, 5.5,
              7.5, 8.5, 10.5, 11.5,
              13.5, 14.5, 16.5, 17.5,
              19.5, 20.5)

for (i in 1:nrow(shal))
{
  arrows(location[i], shal$Mean[i], location[i], shal$Mean[i]+shal$SE[i], code=2, angle=90, length=0.1)
  arrows(location[i], shal$Mean[i], location[i], shal$Mean[i]-shal$SE[i], code=2, angle=90, length=0.1)
} # closes the i in total:1 loop


# Adding statistics

# HC
difference("HC", 15, df3)



text(2, (shal$SE[1] + shal$Mean[1] + 10), "*", cex=2)

# SPON
difference("SPON", 15, df3)

# GORG
difference("GORG", 15, df3)

# MAF
difference("MAF", 15, df3)

# MAC
difference("MAC", 15, df3)

# BC
difference("BC", 15, df3)

# OSUB
difference("OSUB", 15, df3)

# Figure label
legend("topleft", "A - Shallow", bty='n')

##
# Plotting deep
##


barplot(matrix(deep$Mean, nc=length(unique(deep$Category))), beside=TRUE,
        ylim=c(0, 80), axes=FALSE, col=c('grey100', 'grey66'))

axis(1, at=c(2, 5, 8, 11, 14, 17, 20), labels=c("Hard Coral", "Sponge",
                                                "Gorgonian", "Fleshy MA", "Cal. MA", "Black Corals", "Non-Living"),
     cex.axis=0.9, lwd.ticks=FALSE)

axis(2)


# Adding error bars

location <- c(1.5, 2.5, 4.5, 5.5,
              7.5, 8.5, 10.5, 11.5,
              13.5, 14.5, 16.5, 17.5,
              19.5, 20.5)

for (i in 1:nrow(deep))
{
  arrows(location[i], deep$Mean[i], location[i], deep$Mean[i]+deep$SE[i], code=2, angle=90, length=0.1)
  arrows(location[i], deep$Mean[i], location[i], deep$Mean[i]-deep$SE[i], code=2, angle=90, length=0.1)
} # closes the i in total:1 loop


# Adding statistics

# HC
difference("HC", 55, df3)

# SPON
difference("SPON", 55, df3)

# GORG
difference("GORG", 55, df3)
text(8, (deep$SE[5] + deep$Mean[5] + 10), "*", cex=2)


# MAF
difference("MAF", 55, df3)

# MAC
difference("MAC", 55, df3)

# BC
difference("BC", 55, df3)

# OSUB
difference("OSUB", 55, df3)

# Figure label
legend("topleft", "B - MCE", bty='n')
legend(x=16, y=60, c("MPA", "non-MPA"), fill=c('grey100', 'grey66'), cex=0.8)

par(las=0)
mtext("Coverage Type", side=1, line=3.5, outer=TRUE)
mtext("Percentage Cover", side=2, line=0.5, outer=TRUE)


dev.off()






####
# Figure 3 - RDA Patterns - detailed benthic community
####

# set up
rm(list=ls())

library(reshape)
library(vegan)

# Reading in data

df1<- read.csv("ESM 2 Benthic data.csv",
               header=TRUE, colClasses='factor')


##
# Tallying up how many points of each benthic type are recored on each transect
##

# Pasting the site and transect together to create a unqiue
# identifier for each transect

df1$Location <- paste (df1$Site, df1$Depth, df1$Transect, sep="_")

# Combining Cover and Category
df1$Category <- paste(df1$Category, df1$Cover, sep="_")

# Generating frequencies of each location
df2 <- table(c(df1$Location), df1$Category)
df2 <- as.data.frame(df2)
colnames(df2) <- c("Location", "Category", "Frequency")


##
# Calculate percentage cover for each catagory
##

df2$Percentage <- NA

for (i in unique(df2$Location))
{
  
  # Calculating total data points for this location
  temp <- df2[which(df2$Location == i), ]
  total.points <- sum(temp$Frequency)
  
  for (j in unique(df2$Category))
  {
    
    # Calculating the percentage cover of each category
    row.no <- which(df2$Location == i & df2$Category == j)
    df2$Percentage[row.no] <- (df2[row.no, "Frequency"] / total.points) * 100
    
  } # closes the Category loop
} # closes the Location loop


# seperating the location column into site, depth and transect
df2$Site <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 1)
df2$Depth <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 2)
df2$Transect <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 3)


# removing the now redundant location column (first column)
df2 <- subset(df2, select= -1)

##
# Adding in a column with protection status
##

protected <- c("CO", "HE", "PJ", "PT", "SR")
unprotected <- c("PU", "TT", "VB")
df2$Protection <- NA

for (i in 1:nrow(df2))
{
  
  # Protected sites
  if (df2$Site[i] %in% protected)
  {
    df2$Protection[i] <- "PRO"
  }
  
  # Unprotected sites
  if (df2$Site[i] %in% unprotected)
  {
    df2$Protection[i] <- "UNP"
  }
  
} # closes the i in 1:nrow(df3) loop


# Removing non-living benthic habitat types

temp <- as.character(df2$Category)
temp <- strsplit(temp, "_")

temp2 <- c()
for (i in 1:length(temp))
{
  x <- unlist(temp[i])[1]
  temp2 <- c(temp2, x)
}


df2$Category <- temp2

# removing the catagories that are not wanted
df2 <- df2[which(df2$Category!="AINV"), ]
df2 <- df2[which(df2$Category!="OINV"), ]
df2 <- df2[which(df2$Category!="OSUB"), ]

# Casting the dataframe into a format for NMDS
df3 <- cast(df2, Protection + Site + Depth + Transect~ Category, value='Percentage', fun.aggregate=sum)
df3 <- cast(df2, Protection + Site + Depth ~ Category, value='Percentage', fun.aggregate=mean)

# Turning each value into a proportion of the row
df3[, 4:ncol(df3)] <- df3[, 4:ncol(df3)]/rowSums(df3[, 4:ncol(df3)])

# Setting up the figure
tiff("Fig 3 - Benthic RDA.tiff", width = 4.8, height = 9.6, units = 'in', res = 300)
par(mfrow=c(2, 1))
par(mar=c(4, 4, 1, 1))

###
# Shallow reef RDA plot
###

df4 <- df3[which(df3$Depth==15), ]

groups <- df4[, 4:ncol(df4)]
variables <- df4[, 1:2]

# Running redundancy analysis
Shal.NMDS <- rda(groups ~ Protection + Site, data=variables, scale=TRUE)

# Extracting scores
sc <- scores(Shal.NMDS)

# Scaling
y <- 0.5
s1 <- (Shal.NMDS$CCA$eig[1]^y) * sc$sites[, 1]/ (sum(sc$sites[, 1]^2))^y
s2 <- (Shal.NMDS$CCA$eig[2]^y) * sc$sites[, 2]/ (sum(sc$sites[, 2]^2))^y

# Plotting
plot(s1, s2, type='n', xlab=paste("RDA axis 1 (", round(Shal.NMDS$CCA$eig[1]/sum(Shal.NMDS$CCA$eig)*100, 2), "%)", sep=""),
     ylab=paste("RDA axis 2 (", round(Shal.NMDS$CCA$eig[2]/sum(Shal.NMDS$CCA$eig)*100, 2), "%)", sep=""),
     xlim=c(-1, 1), ylim=c(-1, 1))

colour <- "#af8dc3"

character <- c(21, 22)

points(s1, s2, pch=character[as.numeric(as.factor(df4$Protection))],
       bg=colour, col='black', cex=3)


# Sites numbering
# Identifying sites as numbers
sites <- rep(NA, length(df4$Site))

sites[which(df4$Site=='CO')] <- 1
sites[which(df4$Site=='HE')] <- 2
sites[which(df4$Site=='PJ')] <- 3
sites[which(df4$Site=='SR')] <- 4
sites[which(df4$Site=='PT')] <- 5
sites[which(df4$Site=='VB')] <- 6
sites[which(df4$Site=='TT')] <- 7
sites[which(df4$Site=='PU')] <- 8


points(s1, s2, pch=as.character(sites), col='black', cex=1)

# Adding in substrate arrows
arrows(0, 0, sc$species[, 1], sc$species[, 2])

sc$species[4, 2] <- sc$species[4, 2] + 0.1

text(sc$species[3:10, 1], sc$species[3:10, 2], labels=colnames(df4[6:length(df4)]),
     pos=c(1, 2, 1, 4, 3, 1, 3, 1))

# Legend

legend("topleft", "A", bty='n', adj=c(1.5, 0.5), cex=1.5)



###
# Deep reef RDA plot
###

df4 <- df3[which(df3$Depth==55), ]

groups <- df4[, 4:ncol(df4)]
variables <- df4[, 1:2]

# Running redundancy analysis
Deep.NMDS <- rda(groups ~ Protection + Site, data=variables, scale=TRUE)

# Extracting scores
sc <- scores(Deep.NMDS)

# Scaling
y <- 0.5
s1 <- (Deep.NMDS$CCA$eig[1]^y) * sc$sites[, 1]/ (sum(sc$sites[, 1]^2))^y
s2 <- (Deep.NMDS$CCA$eig[2]^y) * sc$sites[, 2]/ (sum(sc$sites[, 2]^2))^y

# Plotting
plot(s1, s2, type='n', xlab=paste("RDA axis 1 (", round(Deep.NMDS$CCA$eig[1]/sum(Deep.NMDS$CCA$eig)*100, 2), "%)", sep=""),
     ylab=paste("RDA axis 2 (", round(Deep.NMDS$CCA$eig[2]/sum(Deep.NMDS$CCA$eig)*100, 2), "%)", sep=""))

colour <- "#7fbf7b"

character <- c(21, 22)

points(s1, s2, pch=character[as.numeric(as.factor(df4$Protection))],
       bg=colour, col='black', cex=3)


# Sites numbering
# Identifying sites as numbers
sites <- rep(NA, length(df4$Site))

sites[which(df4$Site=='CO')] <- 1
sites[which(df4$Site=='HE')] <- 2
sites[which(df4$Site=='PJ')] <- 3
sites[which(df4$Site=='SR')] <- 4
sites[which(df4$Site=='PT')] <- 5
sites[which(df4$Site=='VB')] <- 6
sites[which(df4$Site=='TT')] <- 7
sites[which(df4$Site=='PU')] <- 8


points(s1, s2, pch=as.character(sites), col='black', cex=1)


# Adding in substrate arrows
arrows(0, 0, sc$species[, 1], sc$species[, 2])

text(sc$species[, 1], sc$species[, 2], labels=colnames(df4[4:length(df4)]),
     pos=c(1, 4, 2, 4, 1, 3, 4, 2, 1, 2))

# Legend
legend("topright", legend=c('Shallow', 'MCE', 'MPA', "non-MPA"),
       col=c("#af8dc3", "#7fbf7b", 'black', 'black'), pt.bg='black', pch=c(15, 15, 21, 22),
       ncol=1)


legend("topleft", "B", bty='n', adj=c(1.5, 0.5), cex=1.5)

dev.off()






####
# Figure 4 - Difference in fish richness and biomass with depth inside and outside the MPA
####

# set up
rm(list=ls())

library("vegan")
library("PMCMR")

# Reading in data

df1<- read.csv("ESM 3 Fish biomass.csv",
               header=TRUE)

# Generating unique site and depth combination for each species
df1$combined <- paste(df1$site, df1$depth, df1$transect,
                      df1$family, df1$genus, df1$species,
                      sep="_")

# reducing to one row for each species at each unique location
df2 <- as.data.frame(table(df1$combined))

# returning the location info to original form
df2$site <- unlist(lapply(strsplit(as.character(df2[, 1]), "_"), "[[", 1))
df2$depth <- unlist(lapply(strsplit(as.character(df2[, 1]), "_"), "[[", 2))
df2$transect <- unlist(lapply(strsplit(as.character(df2[, 1]), "_"), "[[", 3))


##
# Total Fish Richness and Abundance per transect
##


# Setting up to loop to count up species
sites <- unique(df2$site)
zones <- sort(as.numeric(unique(df2$depth)))
transects <- sort(as.numeric(unique(df2$transect)))

# results table
df3 <- data.frame(site=NA, zone=NA, transect=NA, richness=NA, abundance=NA, biomass=NA, com.biomass=NA)
# counter for the loop
i <- 1

for (site in sites)
{
  for (zone in zones)
  {
    for (transect in transects)
    {
      # selecting the data
      temp <- df2[which(df2$site==site &
                          df2$depth==zone &
                          df2$transect==transect), ]
      
      # adding in the meta data
      df3[i, c("site", "zone", "transect")] <- c(site, zone, transect)
      
      # richness
      df3[i, "richness"] <- nrow(temp)
      
      # abundance
      df3[i, "abundance"] <- sum(temp$Freq, na.rm=TRUE)
      
      # updating the counter
      i <- i + 1
      
    } # closes the transect in transects loop
  } # closes the zone in zones loop 
}# closes the site in sites loop


##
# Total biomass per transect
##

# loop through the results dataframe

for (i in 1:nrow(df3))
{
  temp <- df1[which(df1$depth==df3$zone[i] &
                      df1$transect==df3$transect[i] &
                      df1$site==df3$site[i]), ]
  
  df3$biomass[i] <- sum(temp$weight.g, na.rm=TRUE)
} # closes the in in 1:nrow(df3)

##
# Total biomass of commercially valuable fish per transect
##


length.weight.table<- read.csv("ESM 4 Commerical value.csv", header=TRUE)
length.weight.table$combined <- paste(length.weight.table$family,
                                      length.weight.table$genus,
                                      length.weight.table$species,
                                      sep="_")

# just commerical value species
df4 <- df1

for(i in 1:nrow(df4))
{
  temp <- as.character(length.weight.table$price[which(length.weight.table$combined == df4$name[i])])
  
  if(length(temp) > 0)
  {
    df4$price[i] <- temp
  } # ends if statement
} # ends for loop


df4 <- df4[which(df4$price=="H" | df4$price=="VH" | df4$price=="M"), ]

# loop through the results dataframe

for (i in 1:nrow(df3))
{
  temp <- df4[which(df4$depth==df3$zone[i] &
                      df4$transect==df3$transect[i] &
                      df4$site==df3$site[i]), ]
  
  df3$com.biomass[i] <- sum(temp$weight.g, na.rm=TRUE)
} # closes the in in 1:nrow(df3)


# Statistics for comparison figure

for (i in 1:nrow(df3))
{
  site <- df3$site[i]
  
  if (site=="TT" | site=="VB" | site=="PU")
  {
    df3$protection[i] <- "OUT"
  } else {
    df3$protection[i] <- "IN"
  }
} # closes the for loop


richness <- aggregate(richness ~ zone + site + protection, data=df3, FUN="mean")

mod1 <- aov(richness ~ zone + protection, data=richness)
plot(mod1)
summary(mod1)
TukeyHSD(mod1)

# overall shallow mean species ricness
mean(richness$richness[which(richness$zone==15)])
# overalll shallow species richness SE
sd(richness$richness[which(richness$zone==15)])/sqrt(8)

# overall MCE mean species ricness
mean(richness$richness[which(richness$zone==55)])
# overalll MCE species richness SE
sd(richness$richness[which(richness$zone==55)])/sqrt(8)


biomass <- aggregate(biomass ~ zone + site + protection, data=df3, FUN="mean")
mod2 <- aov(biomass ~ zone + protection, data=biomass)
summary(mod2)
TukeyHSD(mod2)

mod3 <- aov(biomass ~ protection, data=biomass)

AIC(mod2)
AIC(mod3)

com.biomass <- aggregate(com.biomass ~ zone + site + protection, data=df3, FUN="mean")
mod4 <- aov(com.biomass ~ zone + protection, data=com.biomass)
summary(mod4)
TukeyHSD(mod4)

mod5 <- aov(biomass ~ protection, data=biomass)

AIC(mod4)
AIC(mod5)




###
# Summarising data by site
###

# setting up an empty results table
df4 <- df3[0, ]

sites <- unique(df3$site)
zones <- sort(as.numeric(unique(df3$zone)))
i <- 1

for (site in sites)
{
  for (zone in zones)
  {
    
    # selecting the required data
    temp <- df3[which(df3$site==site & df3$zone==zone),
                c("richness", "abundance", "biomass", "com.biomass")]
    
    # writing the meta-data
    df4[i, c("site", "zone", "transect")] <- c(site, zone, NA)
    
    df4[i, c("richness", "abundance", "biomass", "com.biomass")] <- colMeans(temp)
    
    # update the loop counter
    i <- i + 1
    
  } # closes the zone in zones loop
} # closes the site in sites loop


##
# Sorting the data into inside and outside MPA
##

for (i in 1:nrow(df4))
{
  site <- df4$site[i]
  
  if (site=="TT" | site=="VB" | site=="PU")
  {
    df4$protection[i] <- "OUT"
  } else {
    df4$protection[i] <- "IN"
  }
} # closes the for loop



###
# Calculating the mean richness, abundance and biomass at each depth
###


# results table
df5 <- data.frame(zone=NA, protection=NA, richness=NA, richness.se=NA,
                  abundance=NA, abundance.se=NA,
                  biomass=NA, biomass.se=NA,
                  com.biomass=NA, com.biomass.se=NA)

rownames(df5[1:4, ]) <- 1:4
temp <- expand.grid(zone=unique(df4$zone), protection=unique(df4$protection))

df5$zone <- as.character(temp$zone)
df5$protection <- as.character(temp$protection)
rownames(df5) <- 1:4

# looping through summarising     
for (zone in unique(df5$zone))
{
  for(protection in unique(df5$protection))
  {
    # row placer
    row <- which(df5$zone==zone & df5$protection==protection)
    
    # selecting the correct data
    temp <- df4[which(df4$zone==zone & df4$protection==protection), ]
    
    ###s
    # richness
    df5[row, "richness"] <- sum(temp$richness, na.rm=TRUE)/length(temp$richness)
    
    # richness SE
    df5[row, "richness.se"] <- sqrt(var(temp$richness)/length(temp$richness))
    
    
    ###
    # abundance
    df5[row, "abundance"] <- sum(temp$abundance, na.rm=TRUE)/length(temp$abundance)
    
    # abundance SE
    df5[row, "abundance.se"] <- sqrt(var(temp$abundance)/length(temp$abundance))  
    
    ###
    # biomass
    df5[row, "biomass"] <- sum(temp$biomass, na.rm=TRUE)/length(temp$biomass)
    
    # biomass SE
    df5[row, "biomass.se"] <- sqrt(var(temp$biomass)/length(temp$biomass))    
    ###
    
    # com.biomass
    df5[row, "com.biomass"] <- sum(temp$com.biomass, na.rm=TRUE)/length(temp$com.biomass)
    
    # com.shannon SE
    df5[row, "com.biomass.se"] <- sqrt(var(temp$com.biomass)/length(temp$com.biomass))    
    
  } # closes the protection loop
} # closes the zone loop

###
# Plotting
###

df5 <- df5[c(1, 3, 2, 4), ]

tiff("Fig 4 - Fish richness and biomass.tiff", width = 5, height = 7, units = 'in', res = 300)


par(oma=c(4, 0, 0, 0))
par(mfrow=c(3, 1))
par(mar=c(1, 5, 0.5, 0))
par(cex=1)


# Species Richness

barplot(df5$richness, ylim=c(0, 16), ylab=expression("Richness (150m"^{2} ~")"),
        space=0.25, las=1)

# Adding error bars
total <- nrow(df5)
location <- c(0.75, 2, 3.25, 4.5)
for (i in total:1)
{
  arrows(location[i], df5$richness[i], location[i], df5$richness[i]+df5$richness.se[i], code=2, angle=90, length=0.1)
  arrows(location[i], df5$richness[i], location[i], df5$richness[i]-df5$richness.se[i], code=2, angle=90, length=0.1)
  
} # closes the i in total:1 loop

legend("topleft", "A", bty='n', adj=c(1, -0.5))



# Biomass

barplot(df5$biomass, ylim=c(0, 8000), ylab=expression("Biomass (g/150m"^{2} ~")"),
        space=0.25, las=1)

# Adding error bars
total <- nrow(df5)
location <- c(0.75, 2, 3.25, 4.5)
for (i in total:1)
{
  arrows(location[i], df5$biomass[i], location[i], df5$biomass[i]+df5$biomass.se[i], code=2, angle=90, length=0.1)
  arrows(location[i], df5$biomass[i], location[i], df5$biomass[i]-df5$biomass.se[i], code=2, angle=90, length=0.1)
  
} # closes the i in total:1 loop

legend("topleft", "B", bty='n', adj=c(1, -0.5))


# Commerically Valuable Biomass

barplot(df5$com.biomass, ylim=c(0, 8000), ylab=expression("Biomass (g/150m"^{2} ~")"),
        space=0.25, las=1)

# Adding error bars
total <- nrow(df5)
location <- c(0.75, 2, 3.25, 4.5)
for (i in total:1)
{
  arrows(location[i], df5$com.biomass[i], location[i], df5$com.biomass[i]+df5$com.biomass.se[i], code=2, angle=90, length=0.1)
  arrows(location[i], df5$com.biomass[i], location[i], df5$com.biomass[i]-df5$com.biomass.se[i], code=2, angle=90, length=0.1)
  
} # closes the i in total:1 loop

legend("topleft", "C", bty='n', adj=c(1, -0.5))



axis(1, labels=rep(c("MPA", "non-MPA"), 2), at=location)

mtext("Shallow", side=1, line= 1, outer=TRUE, adj=0.375)
mtext("MCE", side=1, line= 1, outer=TRUE, adj=0.8)

dev.off()





##
# Figure 5 - NMDS of fish community
##

# set up
rm(list=ls())

library('vegan')
library('reshape')

df1<- read.csv("ESM 3 Fish biomass.csv",
               header=TRUE)

# removing species and genus that is NA
df1 <- df1[which(is.na(df1$species)==FALSE), ]
df1 <- df1[which(is.na(df1$genus)==FALSE), ]


# dropping fish with no biomass
df1 <- df1[which(is.na(df1$weight.g)==FALSE), ]



df2 <- cast(df1, site + depth + transect ~ name, value='weight.g', fun.aggregate=sum)
df2 <- as.data.frame(df2)

# turning NaNs into 0
df2[is.na(df2)] <- 0


# Averaging for transects
# back to long form
df3 <- melt(df2, id=c("site", "depth", "transect"))
df3 <- as.data.frame(df3)
# averaging across transects
fish <- cast(df3, site + depth ~ variable, value='value', fun.aggregate=mean)

# Calculating the nMDS
fish.matrix <- vegdist((fish[, 3:ncol(fish)])^(1/4), method='bray')

NMDS <- metaMDS(fish.matrix, k=2)




###
# Plotting
###

# setting the plot layout
tiff("Fig 5 - Fish community NMDS.tiff", width = 6, height = 6, units = 'in', res = 300)


par(mfrow=c(1,1), cex.axis=1.1, cex.lab=1.1)

par(mar=c(4, 4, 1, 1))

# Plot 1 -Adding hard substrate contour lines
plot(NMDS, type='n')

# Identifying protection type
protection <- rep(NA, length(fish$site))

protection[which(fish$site=="PU" | fish$site=="TT" | fish$site=="VB")] <- "OUT"
protection[which(fish$site=="CO" | fish$site=="HE" | fish$site=="PJ" | fish$site=="PT" | fish$site=="SR")] <- "IN"

protection <- c(21, 22)[as.factor(protection)]

# Depths
depths <- c("#af8dc3", "#7fbf7b")[as.factor(fish$depth)]


# Identifying sites as numbers
sites <- rep(NA, length(fish$site))

sites[which(fish$site=='CO')] <- 1
sites[which(fish$site=='HE')] <- 2
sites[which(fish$site=='PJ')] <- 3
sites[which(fish$site=='SR')] <- 4
sites[which(fish$site=='PT')] <- 5
sites[which(fish$site=='VB')] <- 6
sites[which(fish$site=='TT')] <- 7
sites[which(fish$site=='PU')] <- 8


# Adding points
points(NMDS$points[, 1], NMDS$points[, 2],
       pch=protection, bg=depths,  lwd=2, cex=3)

points(NMDS$points[, 1], NMDS$points[, 2], pch=as.character(sites))

# adding the stress label
text(0.35, 0.4, eval(paste("Stress=", round(NMDS$stress, digits=3))), cex=1)

# Legend
legend("bottomright", legend=c('Shallow', 'MCE', 'MPA', "non-MPA"),
       col=c("#af8dc3", "#7fbf7b", 'black', 'black'), pt.bg='black', pch=c(15, 15, 21, 22),
       ncol=1)

dev.off()






####
# Figure 6 - Changes in fish length with depth
####

# set up
#rm(list=ls())

library(sm)
library(KernSmooth)


df1<- read.csv("ESM 5 Fish lengths.csv",
               header=TRUE)


# Removing the point data
df1 <- df1[which((is.na(df1$Length.mm)==FALSE)), ]

# removing the black coral data
df1 <- df1[which(df1$Family!="Antipatharia"), ]

# Preliminary checking
df2 <- df1[, c("Site", "Depth", "Transect", "Length.mm", "Family", "Genus", "Species")]
colnames(df2) <- c("site", "depth", "transect", "length", "family", "genus", "species")

df2$depth <- as.factor(df2$depth)
plot(df2$depth, df2$length)



##
# Compairng shallow fish length distributions between inside and outside the park
##

###
# Shallow
###

# getting data - selecting the sites that are inside and outside
df3 <- df2[which(df2$depth==15), ]

s.in <- df3[which(df3$site=="CO" | df3$site=="SR" | df3$site=="PJ" | df3$site=="PT"| df3$site=="HE"), ]
s.out <- df3[which(df3$site=="TT" | df3$site=="VB" | df3$site=="PU"), ]

# selection of bandwidths
s.in.width <- dpik(s.in$length)
s.out.width <- dpik(s.out$length)

# generating the length distribution outlines
s.in.outline <- bkde(s.in$length, bandwidth=s.in.width)
s.out.outline <- bkde(s.out$length, bandwidth=s.out.width)

# concatinating all lengths
lengths <- c(s.in$length, s.out$length)
groups <- c(rep(1, nrow(s.in)), rep(2, nrow(s.out)))

shal.band <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                nboot=99999, ngrid=500)

shal.n <- length(lengths)


###
# Deep
###

# getting data - selecting the sites that are inside and outside
df3 <- df2[which(df2$depth==55), ]

d.in <- df3[which(df3$site=="CO" | df3$site=="SR" | df3$site=="PJ" | df3$site=="PT"| df3$site=="HE"), ]
d.out <- df3[which(df3$site=="TT" | df3$site=="VB" | df3$site=="PU"), ]

# selection of bandwidths
d.in.width <- dpik(d.in$length)
d.out.width <- dpik(d.out$length)

# generating the length distribution outlines
d.in.outline <- bkde(d.in$length, bandwidth=d.in.width)
d.out.outline <- bkde(d.out$length, bandwidth=d.out.width)

# concatinating all lengths
lengths <- c(d.in$length, d.out$length)
groups <- c(rep(1, nrow(d.in)), rep(2, nrow(d.out)))

deep.band <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                nboot=99999, ngrid=500)

deep.n <- length(lengths)


###
# Running for just the medium and higher priced fish species
###


ref.list <- read.csv("ESM 4 Commerical value.csv", header=TRUE)

# combining family, genus and species into one column
df2$name <- paste(df2$family, df2$genus, df2$species, sep="_")
ref.list$name <- paste(ref.list$family, ref.list$genus, ref.list$species, sep="_")

# Checking the price catagories
unique(ref.list$price)

# Identifying just the higher priced fish species
high.priced <- ref.list[which(ref.list$price=="M" |
                                ref.list$price=="H" |
                                ref.list$price=="VH"), ]

# Selecting just the higher priced fish
high.only <- df2[df2$name %in% high.priced$name, ]

###
# High priced shallow fish
###

# getting data - selecting the sites that are inside and outside
df3 <- high.only[which(high.only$depth==15), ]

h.s.in <- df3[which(df3$site=="CO" | df3$site=="SR" | df3$site=="PJ" | df3$site=="PT"| df3$site=="HE"), ]
h.s.out <- df3[which(df3$site=="TT" | df3$site=="VB" | df3$site=="PU"), ]


# selection of bandwidths
h.s.in.width <- dpik(h.s.in$length)
h.s.out.width <- dpik(h.s.out$length)

# generating the length distribution outlines
h.s.in.outline <- bkde(h.s.in$length, bandwidth=h.s.in.width)
h.s.out.outline <- bkde(h.s.out$length, bandwidth=h.s.out.width)

# concatinating all lengths
lengths <- c(h.s.in$length, h.s.out$length)
groups <- c(rep(1, nrow(h.s.in)), rep(2, nrow(h.s.out)))

h.shal.band <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                  nboot=99999, ngrid=500)

h.shal.n <- length(lengths)


###
# High priced Deep
###

# getting data - selecting the sites that are inside and outside
df3 <- high.only[which(high.only$depth==55), ]

h.d.in <- df3[which(df3$site=="CO" | df3$site=="SR" | df3$site=="PJ" | df3$site=="PT"| df3$site=="HE"), ]
h.d.out <- df3[which(df3$site=="TT" | df3$site=="VB" | df3$site=="PU"), ]


# selection of bandwidths
h.d.in.width <- dpik(h.d.in$length)
h.d.out.width <- dpik(h.d.out$length)

# generating the length distribution outlines
h.d.in.outline <- bkde(h.d.in$length, bandwidth=h.d.in.width)
h.d.out.outline <- bkde(h.d.out$length, bandwidth=h.d.out.width)

# concatinating all lengths
lengths <- c(h.d.in$length, h.d.out$length)
groups <- c(rep(1, nrow(h.d.in)), rep(2, nrow(h.d.out)))

h.deep.band <- sm.density.compare(x=lengths, group=groups, model="equal", method="sj",
                                  nboot=99999, ngrid=500)

h.deep.n <- length(lengths)





###
# Plotting
###


tiff("Fig 6 - Fish length distribution and comparison.tiff", width = 6, height = 6, units = 'in', res = 300)


par(mfrow=c(2, 2))
par(oma=c(3, 3, 0, 0))
par(mar=c(1, 1, 1, 1))

##
# Plotting shallow
##


plot(s.in.outline, lwd=2, type='n', xlim=c(0, 600), ylim=c(0, 0.015),
     axes=FALSE)

axis(1, labels=FALSE)
axis(2, labels=TRUE)

polygon(c(shal.band$eval.points, rev(shal.band$eval.points)), c(shal.band$upper, rev(shal.band$lower)),
        col = "grey", border = NA)

lines(s.in.outline, lwd=2, lty=1)
lines(s.out.outline, lwd=2, lty=2)

text(400, 0.012, eval(paste("n = ", shal.n)))
text(400, 0.010, "p < 0.001")
text(0, 0.014, "A - All shallow reef fishes", pos=4, cex=0.8)


##
# Plotting Deep
##


plot(d.in.outline, lwd=2, type='n', xlim=c(0, 600), ylim=c(0, 0.015),
     axes=FALSE)


axis(1, labels=FALSE)
axis(2, labels=FALSE)

polygon(c(deep.band$eval.points, rev(deep.band$eval.points)), c(deep.band$upper, rev(deep.band$lower)),
        col = "grey", border = NA)

lines(d.in.outline, lwd=2, lty=1)
lines(d.out.outline, lwd=2, lty=2)

text(400, 0.012, eval(paste("n = ", deep.n)))
text(400, 0.010, eval(paste("p = ", round(deep.band$p, 3))))
text(0, 0.014, "B - All MCE fishes", pos=4, cex=0.8)

legend(200, 0.009, legend=c("MPA", "non-MPA"), lty=c(1, 2), lwd=2, bty="n")



####
# Plotting high value shallow
####


plot(h.s.in.outline, lwd=2, type='n', xlim=c(0, 600), ylim=c(0, 0.015),
     axes=FALSE)


axis(1, labels=TRUE)
axis(2, labels=TRUE)

polygon(c(h.shal.band$eval.points, rev(h.shal.band$eval.points)), c(h.shal.band$upper, rev(h.shal.band$lower)),
        col = "grey", border = NA)

lines(h.s.in.outline, lwd=2, lty=1)
lines(h.s.out.outline, lwd=2, lty=2)

text(400, 0.012, eval(paste("n = ", h.shal.n)))
text(400, 0.010, "p < 0.001")
text(0, 0.014, "C - Commercially valuable shallow reef fishes", pos=4, cex=0.8)




##
# Plotting High value deep
##


plot(h.d.in.outline, lwd=2, type='n', xlim=c(0, 600), ylim=c(0, 0.015),
     axes=FALSE)


axis(1, labels=TRUE)
axis(2, labels=FALSE)

polygon(c(h.deep.band$eval.points, rev(h.deep.band$eval.points)), c(h.deep.band$upper, rev(h.deep.band$lower)),
        col = "grey", border = NA)

lines(h.d.in.outline, lwd=2, lty=1)
lines(h.d.out.outline, lwd=2, lty=2)

text(400, 0.012, eval(paste("n = ", h.deep.n)))
text(400, 0.010, eval(paste("p = ", round(h.deep.band$p, 3))))
text(0, 0.014, "D - Commercially valuable MCE fishes", pos=4, cex=0.8)



# Axis labels

mtext("Fish length (mm)", 1, line=1.5, outer=TRUE)
mtext("Probability Density", 2, line=1.5, outer=TRUE)


dev.off()






####
# Table 1 - Benthic community PERMANOVA
####

# set up
rm(list=ls())

library(reshape)
library(vegan)

# Reading in data

df1<- read.csv("ESM 2 Benthic data.csv",
               header=TRUE, colClasses='factor')


##
# Tallying up how many points of each benthic type are recored on each transect
##

# Pasting the site and transect together to create a unqiue
# identifier for each transect

df1$Location <- paste (df1$Site, df1$Depth, df1$Transect, sep="_")

# Combining Cover and Category
df1$Category <- paste(df1$Category, df1$Cover, sep="_")

# Generating frequencies of each location
df2 <- table(c(df1$Location), df1$Category)
df2 <- as.data.frame(df2)
colnames(df2) <- c("Location", "Category", "Frequency")


##
# Calculate percentage cover for each catagory
##

df2$Percentage <- NA

for (i in unique(df2$Location))
{
  
  # Calculating total data points for this location
  temp <- df2[which(df2$Location == i), ]
  total.points <- sum(temp$Frequency)
  
  for (j in unique(df2$Category))
  {
    
    # Calculating the percentage cover of each category
    row.no <- which(df2$Location == i & df2$Category == j)
    df2$Percentage[row.no] <- (df2[row.no, "Frequency"] / total.points) * 100
    
  } # closes the Category loop
} # closes the Location loop


# seperating the location column into site, depth and transect
df2$Site <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 1)
df2$Depth <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 2)
df2$Transect <- sapply(strsplit(as.character(df2$Location), '_'), "[[", 3)


# removing the now redundant location column (first column)
df2 <- subset(df2, select= -1)

##
# Adding in a column with protection status
##

protected <- c("CO", "HE", "PJ", "PT", "SR")
unprotected <- c("PU", "TT", "VB")
df2$Protection <- NA

for (i in 1:nrow(df2))
  {
  
  # Protected sites
  if (df2$Site[i] %in% protected)
    {
    df2$Protection[i] <- "PRO"
    }
  
  # Unprotected sites
  if (df2$Site[i] %in% unprotected)
    {
    df2$Protection[i] <- "UNP"
    }
  
  } # closes the i in 1:nrow(df3) loop

# Removing non-living benthic habitat types and simplifying down Categorys
temp <- as.character(df2$Category)
temp <- strsplit(temp, "_")

temp2 <- c()
  for (i in 1:length(temp))
  {
  x <- unlist(temp[i])[1]
  temp2 <- c(temp2, x)
  }


df2$Category <- temp2

# removing the catagories that are not wanted
df2 <- df2[which(df2$Category!="AINV"), ]
df2 <- df2[which(df2$Category!="OINV"), ]
df2 <- df2[which(df2$Category!="OSUB"), ]

  

# Casting the dataframe into a format for PERMANOVA
df3 <- cast(df2, Protection + Site + Depth + Transect ~ Category, value='Percentage', fun.aggregate=sum)

# Adding lattitude of survey sites
gps.points <- read.csv("ESM 1 Cozumel GPS points.csv")

for (i in 1:nrow(df3))
{
  row.no <- which(gps.points$Site==df3$Site[i] & gps.points$Depth==df3$Depth[i])
  df3$lattitude[i] <- gps.points$Latitude[row.no]
} # closes the for loop



benthic.matrix <- vegdist(df3[, 5:(ncol(df3)-1)], method="bray")

# Running PERMANOVA

model <- adonis(benthic.matrix ~ df3$lattitude + df3$Protection * df3$Site * df3$Depth,
                distance="bray", permutations = 99999)

results <- as.data.frame(model$aov.tab)



write.csv(results, file="Table_benthic_permanova.csv", row.names=TRUE)






##
# Table 2 - Fish community PERMANOVA
##

# set up
rm(list=ls())

library('vegan')
library('reshape')

# Reading in data

df1<- read.csv("ESM 3 Fish biomass.csv",
               header=TRUE)

# removing species and genus that is NA
df1 <- df1[which(is.na(df1$species)==FALSE), ]
df1 <- df1[which(is.na(df1$genus)==FALSE), ]


# dropping fish with no biomass
df1 <- df1[which(is.na(df1$weight.g)==FALSE), ]



df2 <- cast(df1, site + depth + transect ~ name, value='weight.g', fun.aggregate=sum)
df2 <- as.data.frame(df2)

# turning NaNs into 0
df2[is.na(df2)] <- 0


fish <- df2

# Calculating the matrix
fish.matrix <- vegdist((fish[, 4:ncol(fish)])^(1/4), method='bray')

# Allocating the protection status
protection <- rep(NA, length(fish$site))

protection[which(fish$site=="PU" | fish$site=="TT" | fish$site=="VB")] <- "OUT"
protection[which(fish$site=="CO" | fish$site=="HE" | fish$site=="PJ" | fish$site=="PT" | fish$site=="SR")] <- "IN"

# Adding lattitude of survey sites
gps.points <- read.csv("ESM 1 Cozumel GPS points.csv")

latitude <- rep(NA, length(fish$site))

for (i in 1:length(latitude))
  {
  latitude[i] <- gps.points$Latitude[which(gps.points$Site==fish$site[i] & gps.points$Depth==fish$depth[i])]
  } # closes the for loop

# Fitting the model
model <- adonis(fish.matrix ~ latitude + protection * fish$site * fish$depth,
                distance="bray", permutations = 99999)

results <- as.data.frame(model$aov.tab)

write.csv(results, file="Table_fish_permanova.csv", row.names=TRUE)





####
# Table 3 - Fish species correlation with PCOA
###

# set up
rm(list=ls())

library('vegan')
library('reshape')
library('dplyr')

# Reading in data

df1<- read.csv("ESM 3 Fish biomass.csv",
               header=TRUE)

# removing species and genus that is NA
df1 <- df1[which(is.na(df1$species)==FALSE), ]
df1 <- df1[which(is.na(df1$genus)==FALSE), ]


# Calculating biomass


# dropping fish with no biomass
df1 <- df1[which(is.na(df1$weight.g)==FALSE), ]



df2 <- cast(df1, site + depth + transect ~ name, value='weight.g', fun.aggregate=sum)
df2 <- as.data.frame(df2)

# turning NaNs into 0
df2[is.na(df2)] <- 0


# Averaging for transects
# back to long form
df3 <- melt(df2, id=c("site", "depth", "transect"))
df3 <- as.data.frame(df3)
# averaging across transects
fish <- cast(df3, site + depth ~ variable, value='value', fun.aggregate=mean)


# Protection status
protection <- rep(NA, length(fish$site))

protection[which(fish$site=="PU" | fish$site=="TT" | fish$site=="VB")] <- "OUT"
protection[which(fish$site=="CO" | fish$site=="HE" | fish$site=="PJ" | fish$site=="PT" | fish$site=="SR")] <- "IN"


# Running PCOA
fish.matrix <- vegdist((fish[, 3:ncol(fish)])^(1/4), method='bray')

PCOA <- cmdscale(fish.matrix)


# Extracting PCOA scores
x <- scores(PCOA)[, 1]
y <- scores(PCOA)[, 2]



## Testing correlations with PCOA axis

df5 <- aggregate(value ~ variable + site + depth, data=df3,FUN=mean)


# Blank vector to store taxonmic groups in that correlate with PCOA
species <- c()

# correlation threshold value
threashold <- 0.4

# First axis
for (i in unique(df5$variable))
  {
  temp <- df5[which(df5$variable==i), ]
  
  temp <- left_join(fish[, 1:2], temp)
  
  cor.result <- cor.test(temp$value, x, method="pearson")
  if (abs(cor.result$estimate) > threashold)
    {
    species <- c(species, i)
    } # closes if
  
  }

x.species <- species


# Second axis
species <- c()

for (i in unique(df5$variable))
  {
  temp <- df5[which(df5$variable==i), ]
  
  temp <- left_join(fish[, 1:2], temp)
  
  cor.result <- cor.test(temp$value, y, method="pearson")
  if (abs(cor.result$estimate) > threashold)
    {
    species <- c(species, i)
    } # closes if
  
  }
y.species <- species


##
# Recording differences in abundance between depth
##

# Identifying protection type for the site level data
protection <- rep(NA, length(df5$site))

protection[which(df5$site=="PU" | df5$site=="TT" | df5$site=="VB")] <- "OUT"
protection[which(df5$site=="CO" | df5$site=="HE" | df5$site=="PJ" | df5$site=="PT" | df5$site=="SR")] <- "IN"

df5$protection <- protection

# Identifying protection type for the transect level data
protection <- rep(NA, length(df3$site))

protection[which(df3$site=="PU" | df3$site=="TT" | df3$site=="VB")] <- "OUT"
protection[which(df3$site=="CO" | df3$site=="HE" | df3$site=="PJ" | df3$site=="PT" | df3$site=="SR")] <- "IN"

df3$protection <- protection


# Adding lattitude of survey sites
gps.points <- read.csv("ESM 1 Cozumel GPS points.csv")

latitude <- rep(NA, length(df3$site))

for (i in 1:length(latitude))
{
  latitude[i] <- gps.points$Latitude[which(gps.points$Site==df3$site[i] & gps.points$Depth==df3$depth[i])]
} # closes the for loop

df3$latitude <- latitude

# species list to store correlating species

species.list <- c(x.species, y.species)

# results dataframe
results.x <- data.frame(group=species.list,
                        in.15=NA, in.15.se=NA,
                        out.15=NA, out.15.se=NA,
                        in.55=NA, in.55.se=NA,
                        out.55=NA, out.55.se=NA,
                        f.d=NA, p.d=NA,
                        f.p=NA, p.p=NA,
                        f.i=NA, p.i=NA)

for (i in 1:length(species.list))
{
  species <- species.list[i]
  
  results.x[i, "in.15"] <- mean(df5[which(df5$depth=="15" & df5$protection=="IN" & df5$variable == species), 'value'])
  results.x[i, "in.15.se"] <- sd(df5[which(df5$depth=="15" & df5$protection=="IN" & df5$variable == species), 'value']) / sqrt(length(df5[which(df5$depth=="15" & df5$protection=="IN" & df5$variable == species), 'value']))
  
  results.x[i, "out.15"] <- mean(df5[which(df5$depth=="15" & df5$protection=="OUT" & df5$variable == species), 'value'])
  results.x[i, "out.15.se"] <- sd(df5[which(df5$depth=="15" & df5$protection=="OUT" & df5$variable == species), 'value']) / sqrt(length(df5[which(df5$depth=="15" & df5$protection=="OUT" & df5$variable == species), 'value']))

  results.x[i, "in.55"] <- mean(df5[which(df5$depth=="55" & df5$protection=="IN" & df5$variable == species), 'value'])
  results.x[i, "in.55.se"] <- sd(df5[which(df5$depth=="55" & df5$protection=="IN" & df5$variable == species), 'value']) / sqrt(length(df5[which(df5$depth=="55" & df5$protection=="IN" & df5$variable == species), 'value']))
  
  results.x[i, "out.55"] <- mean(df5[which(df5$depth=="55" & df5$protection=="OUT" & df5$variable == species), 'value'])
  results.x[i, "out.55.se"] <- sd(df5[which(df5$depth=="55" & df5$protection=="OUT" & df5$variable == species), 'value']) / sqrt(length(df5[which(df5$depth=="55" & df5$protection=="OUT" & df5$variable == species), 'value']))
  

  temp <- df3[which(df3$variable == species), ]
  
  temp <- adonis(vegdist(temp$value, method="euclid") ~ temp$latitude + temp$protection * temp$site * temp$depth,
                 method="euclid", permutations = 99999)
  
  
  
  #depth
  results.x[i, "f.d"] <- temp$aov.tab[4,4]
  results.x[i, "p.d"] <- temp$aov.tab[4,6]
  
  #protection
  results.x[i, "f.p"] <- temp$aov.tab[2,4]
  results.x[i, "p.p"] <- temp$aov.tab[2,6]
  
  #interaction
  results.x[i, "f.i"] <- temp$aov.tab[5,4]
  results.x[i, "p.i"] <- temp$aov.tab[5,6]
} # closes the results for loop


# writing results to a table
write.csv(results.x, "Table_3_fish_species_interaction.csv", row.names=FALSE)