rm(list=ls())
if (!require("scales")) install.packages("scales"); library("scales")
if (!require("vegan")) install.packages("vegan"); library("vegan")
if (!require("tidyverse")) install.packages("tidyverse"); library("tidyverse")
if (!require("gridExtra")) install.packages("gridExtra"); library("gridExtra")
palette <- hue_pal()(11)

fish_multivariate <- read.csv("fish_multivariate.csv")
inverts_multivariate <- read.csv("inverts_multivariate.csv")
benthic_multivariate <- read.csv("benthic_multivariate.csv")

# dbRDA to predict multivariate space with PC1 and PC2, shore category and protection level for fish communities
fishspecies <- decostand(fish_multivariate[,8:69], method = "total")
environment <- fish_multivariate[,3:6]
dbRDA_fish <- capscale(fishspecies~PC1+PC2+Category+Protection, environment, dist = "bray")
anova(dbRDA_fish, by = "terms", permu = 1000) # test variables for sig.
anova(dbRDA_fish, by="axis", perm.max=1000) # test axes for significance
anova(dbRDA_fish)
RsquareAdj(dbRDA_fish)
# extract % explained by first two CAP axes
perc <- round(100*(summary(dbRDA_fish)$cont$importance[2, 1:2]), 2)
## extract scores - these are coordinates in the RDA space
sc_si <- scores(dbRDA_fish, display="sites", choices=c(1,2), scaling=2)
sc_sp <- scores(dbRDA_fish, display="species", choices=c(1,2), scaling=2)
sc_bp <- scores(dbRDA_fish, display="bp", choices=c(1, 2), scaling=2)

#get top 5 most influential species
fish_sp <- sc_sp
# vector lengths
fish_vec_len <- sqrt(fish_sp[,1]^2 + fish_sp[,2]^2)
# pick top 5 species
fishtop5 <- names(sort(fish_vec_len, decreasing = TRUE))[1:5]
fish_sp_top5 <- fish_sp[fishtop5, ]

# triplot based on the ordiplot() format
par(mar = c(5.1,5.1,4.1,2.1))
fishbiplot <- ordiplot(dbRDA_fish, scaling = 2, type = "none", cex.lab = 1.5, cex.axis = 1.25, xlim = c(-4, 4), ylim = c(-5,5), xlab = paste0("CAP1 (", perc[1], "%)"), ylab = paste0("CAP2 (", perc[2], "%)"), xpd = TRUE)
points(fishbiplot, "sites", pch = c(21,24)[as.factor(fish_multivariate$Category)], cex = 1.5, col = "black", bg = alpha(palette[as.factor(fish_multivariate$Location)],0.8))
#text(b, "species", col = "darkslategrey", cex = 1.5)
arrows(0,0, # start them from (0,0)
       sc_bp[,1]*3, sc_bp[,2]*3, # end them at the score value
       col = "black")
text(x = sc_bp[,1]*3.25, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2]*3.25,
     labels = rownames(sc_bp),
     col = "black",
     cex = 1.5,
     font = 2)

text(fish_sp_top5[,1]*1.5, fish_sp_top5[,2]*1.5,
     labels = rownames(fish_sp_top5),
     col = "darkslategrey",
     cex = 1.5)
legend("bottomleft", unique(sort(fish_multivariate$Location)), pt.bg = palette[unique(as.factor(fish_multivariate$Location))], pch = c(24,21,21,24,24,21,21,21,24,24,24), bty = "n", cex =1.5, xpd = TRUE)

fishspecies_scores <- round(scores(dbRDA_fish)$species,3)


# dbRDA to predict multivariate space with PC1 and PC2, shore category and protection level for inverts communities

### inverts -----------------

# dbRDA to predict multivariate space with PC1 and PC2--------------
invertspecies <- decostand(inverts_multivariate[8:36], method = "total")
environment <- inverts_multivariate[,c(3:6)]
invertspecies0001 <- invertspecies +0.0001 # adding a very small number to the 0
dbRDA_inverts <- capscale(invertspecies0001~PC1+PC2+Category+Protection, environment, dist = "bray")
anova(dbRDA_inverts, by = "terms", permu = 1000) # test variables for sig.
anova(dbRDA_inverts, by="axis", perm.max=1000) # test axes for significance
anova(dbRDA_inverts)
RsquareAdj(dbRDA_inverts)
# extract % explained by first two CAP axes
perc <- round(100*(summary(dbRDA_inverts)$cont$importance[2, 1:2]), 2)
## extract scores - these are coordinates in the RDA space
sc_si <- scores(dbRDA_inverts, display="sites", choices=c(1,2), scaling=2)
sc_sp <- scores(dbRDA_inverts, display="species", choices=c(1,2), scaling=2)
sc_bp <- scores(dbRDA_inverts, display="bp", choices=c(1, 2), scaling=2)

#get top 5 most influential species
inverts_sp <- sc_sp
# vector lengths
inverts_vec_len <- sqrt(inverts_sp[,1]^2 + inverts_sp[,2]^2)
# pick top 5 species
invertstop5 <- names(sort(inverts_vec_len, decreasing = TRUE))[1:5]
inverts_sp_top5 <- inverts_sp[invertstop5, ]


# triplot based on the ordiplot() format
par(mar = c(5.1,5.1,4.1,2.1))
invertsbiplot <- ordiplot(dbRDA_inverts, scaling = 2, type = "none", cex.lab = 1.5, cex.axis = 1.25, xlim = c(-5, 5), ylim = c(-5,5), xlab = paste0("CAP1 (", perc[1], "%)"), ylab = paste0("CAP2 (", perc[2], "%)"), xpd = TRUE)
points(invertsbiplot, "sites", pch = c(21,24)[as.factor(inverts_multivariate$Category)], cex = 1.5, col = "black", bg = alpha(palette[as.factor(inverts_multivariate$Location)],0.8))
#text(b, "species", col = "darkslategrey", cex = 1.5)
arrows(0,0, # start them from (0,0)
       sc_bp[,1]*3, sc_bp[,2]*3, # end them at the score value
       col = "black")
text(x = sc_bp[,1]*3.25, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2]*3.25, 
     labels = rownames(sc_bp), 
     col = "black", 
     cex = 1.5,
     font = 2)
#legend("bottomleft", unique(sort(inverts_multivariate$Location)), pt.bg = palette[unique(as.factor(inverts_multivariate$Location))], pch = c(24,21,21,24,24,21,21,21,24,24,24), bty = "n", cex =1, xpd = TRUE)

text(inverts_sp_top5[,1]*1.25, inverts_sp_top5[,2]*1.25,
     labels = rownames(inverts_sp_top5),
     col = "darkslategrey",
     cex = 1.5)

invertspecies_scores <- round(scores(dbRDA_inverts)$species,3)


# dbRDA to predict multivariate space with PC1 and PC2, shore category and protection level for the benthic communities
species <- benthic_multivariate[,c(10:25)]
environment <- benthic_multivariate[,c(4:7)]
# running dbRDA
dbRDA <- capscale(species~PC1+PC2+Category+Protection, environment, dist = "bray")
anova(dbRDA, by = "terms", permu = 1000) # test variables for sig.
anova(dbRDA, by="axis", perm.max=1000) # test axes for significance
anova(dbRDA)
RsquareAdj(dbRDA)
# extract % explained by first two CAP axes
perc <- round(100*(summary(dbRDA)$cont$importance[2, 1:2]), 2)
## extract scores - these are coordinates in the RDA space
sc_si <- scores(dbRDA, display="sites", choices=c(1,2), scaling=2)
sc_sp <- scores(dbRDA, display="species", choices=c(1,2), scaling=2)
sc_bp <- scores(dbRDA, display="bp", choices=c(1, 2), scaling=2)

#get top 5 most influential species
coral_sp <- sc_sp
# vector lengths
coral_vec_len <- sqrt(coral_sp[,1]^2 + coral_sp[,2]^2)
# pick top 5 species
coraltop5 <- names(sort(coral_vec_len, decreasing = TRUE))[1:5]
coral_sp_top5 <- coral_sp[coraltop5, ]

# triplot based on the ordiplot() format
par(mar = c(5.1,5.1,4.1,2.1))
a <- ordiplot(dbRDA, scaling = 2, type = "none", cex.lab = 1.5, cex.axis = 1.25, xlim = c(-5, 5), ylim = c(-5,5), xlab = paste0("CAP1 (", perc[1], "%)"), ylab = paste0("CAP2 (", perc[2], "%)"), xpd = TRUE)
points(a, "sites", pch = c(21,24)[as.factor(benthic_multivariate$Category)], cex = 1.5, col = "black", bg = alpha(palette[as.factor(benthic_multivariate$Location)],0.8))
#text(a, "species", col = "darkslategrey", cex = 1.5) #######ISSUE check the column headers are good 
arrows(0,0, # start them from (0,0)
       sc_bp[,1]*3, sc_bp[,2]*3, # end them at the score value
       col = "black")
text(x = sc_bp[,1]*3.25, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2]*3.25, 
     labels = rownames(sc_bp), 
     col = "black", 
     cex = 1.5,
     font = 2)
text(coral_sp_top5[,1]*1.5, coral_sp_top5[,2]*1.5,
     labels = rownames(coral_sp_top5),
     col = "darkslategrey",
     cex = 1.5)


legend("bottomleft", unique(sort(benthic_multivariate$Location)), pt.bg = palette[unique(as.factor(benthic_multivariate$Location))], pch = c(24,21,21,24,24,21,21,21,24,24,24), bty = "n", cex =1.25, xpd = TRUE)
benthicspecies_scores <- round(scores(dbRDA)$species,3)