library(readr)
library(vegan)
library(tidyverse)
library(usdm)
library(adespatial)
library(adegraphics)
setwd(".../Naas_T")


data <- read_csv("Naas_T-data.csv",
                 col_types = cols(`Wetland type` =
                                    col_factor(levels = c("Temp vlei",
                                                        "Large dam", "Small dam", "River edge",
                                                        "Fynbos pool")), temporary = col_factor(levels = c("temp",
                                                                                                           "perm")), Fish = col_factor(levels = c("no-fish",
                                                                                                                                                  "fish"))))
data #wetland type, hydroperiod, and fish treated as factors

coords <- data[,4:5]
sp_compos <- data %>%
  dplyr::select(Xl:Hh)
env <- data %>%
  dplyr::select(c(`Wetland type`, temporary, area, perimeter, pH, Conductivity,Fish ))
env



# presence-absence of frog species only
pa <- decostand(sp_compos, "pa")


plot(data[,4:5])


mantel.correlog(vegdist(sp_compos), XY = coords)#no spatial autocorrelation
plot(mantel.correlog(vegdist(sp_compos), XY = coords))
mems <- dbmem(data[,4:5], MEM.autocor = "positive")
#building MEMs and only taking eigenvectors with positive Moran's I
s.value(coords, mems[,1:2], include.origin = FALSE)

MEMs <- cbind(MEM1=mems$MEM1, MEM2=mems$MEM2)


gm.pa <- rda(pa~`Wetland type`+ temporary+ area+ perimeter+ pH+ Conductivity+Fish + Condition(MEMs), data=env)
gm.pa

## “The response variables explain 27.4% of the variation in anuran community composition across sites, while site control variables explains 9.0% of the variation.”
RsquareAdj(gm.pa)$adj.r.squared
anova.cca(gm.pa, step = 10000)


## The model explains 10.0% of the variation in frog community composition across sites. It is also statistically significant (p = 0.002)


# A description of the significance of environmental variables:

varpart(sp_compos, env, MEMs)

anova.cca(gm.pa, by = "term")

anova.cca(gm.pa, by = "axis")

RsquareAdj(gm.pa)

#  Marginal effects of terms for the environmental variables from the partial RDA model.

anova.cca(gm.pa, step = 10000, by = "margin")

# sig differences are driven by fish - no fish & fynbos pools and small dams. 
anova.cca(gm.pa, step = 10000, by = "onedf")




#######################################################
####################### Figure 2 ######################
#######################################################

# Set up a blank plot with scaling, axes, and labels
plot(gm.pa,
     scaling = 2, # set scaling type 
     type = "none", # this excludes the plotting of any points from the results
     frame = FALSE,
     # set axis limits
     xlim = c(-1,1), 
     ylim = c(-1,1),
     # label the plot (title, and axes)
     main = "Triplot Frog RDA - scaling 2",
     xlab = paste0("RDA1 (", perc[1], "%)"), 
     ylab = paste0("RDA2 (", perc[2], "%)") 
)
perc <- round(100*(summary(gm.pa)$cont$importance[2, 1:2]), 2)
sc_si <- scores(gm.pa, display="sites", choices=c(1,2), scaling=2)
sc_sp <- scores(gm.pa, display="species", choices=c(1,2), scaling=2)
sc_bp <- scores(gm.pa, display="bp", choices=c(1, 2), scaling=2)
sc_sp

# add points for site scores
points(sc_si, 
       pch = 21, # set shape (here, circle with a fill colour)
       col = "black", # outline colour
       bg = "steelblue", # fill colour
       cex = 1.2) # size
# add points for species scores
points(sc_sp, 
       pch = 22, # set shape (here, square with a fill colour)
       col = "black",
       bg = "#f2bd33", 
       cex = 1.2)
# add text labels for species abbreviations
text(sc_sp + c(0.03, 0.09), # adjust text coordinates to avoid overlap with points 
     labels = rownames(sc_sp), 
     col = "grey40", 
     font = 2, # bold
     cex = 0.6)
# add arrows for effects of the explanatory variables
arrows(0,0, # start them from (0,0)
       sc_bp[,1], sc_bp[,2], # end them at the score value
       col = "red", 
       lwd = 2)
# add text labels for arrows
text(x = sc_bp[,1] -0.1, # adjust text coordinate to avoid overlap with arrow tip
     y = sc_bp[,2] - 0.03, 
     labels = rownames(sc_bp), 
     col = "red", 
     cex = 1, 
     font = 1.5)





############################ Analyses on reduced dataset of permanent water bodies ###############################
# Because invasive fish can only occur in permanent water bodies, we carried out an adonis analysis on fish within permanent water only.
##################################################################################################################

data <- read.csv("Naas_T-data.csv")
frog_data <- data[,19:29]
va_data <- data[,1:17]
var_data <- va_data[, -c(6, 9, 10, 11, 13, 14)]


data_nf <- subset(data, var_data$temporary=="perm")
frog_data_nf <- data_nf[,19:29] 
va_data_nf <- data_nf[,4:17]

summary(frog_data_nf)
summary(va_data_nf)

var_data_nf <- va_data_nf[, -c(6, 9, 10, 11)]
frog_pa_nf <- decostand(frog_data_nf, "pa")


# We used Jaccard dissimilarity because of the presence absence nature of data. We tried k2, 3 and 4
pa.nf.nmds.k2 <- metaMDS(frog_pa_nf, trymax = 100, k = 2, distance = "jaccard", autotransform = F)
pa.nf.nmds.k3 <- metaMDS(frog_pa_nf, trymax = 100, k = 3, distance = "jaccard", autotransform = F)
pa.nf.nmds.k4 <- metaMDS(frog_pa_nf, trymax = 100, k = 4, distance = "jaccard", autotransform = F)

#Test to see how many dimensions the nmds should have
pa.nf.nmds.k2$stress
pa.nf.nmds.k3$stress
pa.nf.nmds.k4$stress


## Shepard Plot method for class 'metaMDS'
stressplot(pa.nf.nmds.k2)
stressplot(pa.nf.nmds.k3)
stressplot(pa.nf.nmds.k4)



## Default S3 method:
stressplot(object, dis, pch, p.col = "blue", l.col = "red", 
           lwd = 2, ...) 

nfismorethan0.05 <- (pa.nf.nmds.k4$stress - pa.nf.nmds.k3$stress)
nfismorethan0.05
#  2D is better than 3D as stress drops by 0.043. Critical change acceptable is 0.05


rm(pa.nf.nmds.k3)
rm(pa.nf.nmds.k4)


fish.perm.nf <- adonis2(frog_pa_nf ~ var_data_nf$Fish, data = var_data_nf)  
fish.perm.nf
coef(fish.perm.nf)
summary(fish.perm.nf)
print(fish.perm.nf)



fort.nf <- fortify(pa.nf.nmds.k2)


NMDS_nf <- ggplot(pa.nf.nmds.k2, aes(x = NMDS1, y = NMDS2)) + geom_point(aes(NMDS1, NMDS2, colour = var_data_nf$Fish), size = 4) + labs(colour = "Presence or absence of invasive fish")
NMDS_nf

###################### Figure 3 #############################
# plot NMDS of sites with permanent water with or without invasive fish
#############################################################

scn_si <- scores(pa.nf.nmds.k2, display="sites", choices=c(1,2), scaling=2)
scn_sp <- scores(pa.nf.nmds.k2, display="species", choices=c(1,2), scaling=2)
scn_bp <- scores(pa.nf.nmds.k2, display="species", choices=c(1,2), scaling=2)

scores(pa.nf.nmds.k2)
pa.nf.nmds.k2

# Set up a blank plot with scaling, axes, and labels
plot(pa.nf.nmds.k2,
     scaling = 2, # set scaling type 
     type = "none", # this excludes the plotting of any points from the results
     frame = FALSE,
     # set axis limits
     xlim = c(-1,1), 
     ylim = c(-1,1),
     # label the plot (title, and axes)
     main = "Triplot Frog NMDS - scaling 2",
     xlab = paste0("NMDS1"), 
     ylab = paste0("NMDS2") 
)

# add points for site scores
points(scn_si, 
       pch = 21, # set shape (here, circle with a fill colour)
       col = "black", # outline colour
       bg = "steelblue", # fill colour
       cex = 1.2) # size
# add points for species scores
points(scn_sp, 
       pch = 22, # set shape (here, square with a fill colour)
       col = "black",
       bg = "#f2bd33", 
       cex = 1.2)
# add text labels for species abbreviations
text(scn_sp + c(0.03, 0.09), # adjust text coordinates to avoid overlap with points 
     labels = rownames(sc_sp), 
     col = "grey40", 
     font = 2, # bold
     cex = 0.6)
# add arrows for effects of the explanatory variables
arrows(0,0, # start them from (0,0)
       scn_bp[,1], scn_bp[,2], # end them at the score value
       col = "red", 
       lwd = 2)
# add text labels for arrows
text(x = scn_bp[,1] -0.1, # adjust text coordinate to avoid overlap with arrow tip
     y = scn_bp[,2] - 0.03, 
     labels = rownames(scn_bp), 
     col = "red", 
     cex = 1, 
     font = 1.5)


###################### Figure 3 #############################
# plot NMDS of sites with permanent water with or without invasive fish

Fig3 <- ggplot() +
  stat_ellipse(data = subset(fort.nf, Score == 'sites'), 
               aes(x = NMDS1, y = NMDS2, 
                   fill=var_data_nf$Fish,
                   colour = var_data_nf$Fish), 
               geom="polygon", level = 0.5, alpha=0.2) +
  geom_point(data = subset(fort.nf, Score == 'sites'),
             mapping = aes(x = NMDS1, y = NMDS2, colour = var_data_nf$Fish),
             alpha = 0.5)+
  geom_segment(data = subset(fort.nf, Score == 'species'),
               mapping = aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
               arrow = arrow(length = unit(0.015, "npc"),
                             type = "closed"),
               colour = "darkgray",
               size = 0.4)+
  geom_text(data = subset(fort.nf,Score == 'species'),
            mapping = aes(label = Label, x = NMDS1 * 1.1, y = NMDS2 * 1.1)) +
  geom_abline(intercept = 0,slope = 0,linetype = "dashed", size=0.8, colour="gray") +
  geom_vline(aes(xintercept=0), linetype ="dashed", size=0, colour="gray") +
  theme(legend.position = "none")+
  coord_fixed(xlim = c(-1,1), ylim = c(-1,1))

Fig3

ggsave("Fig3.pdf", width = 75, height = 75, units = "mm")
ggarrange(Fig3, ncol = 1,nrow = 1)
dev.off()


########################################################################################################################
###### We acknowledge the help of Diogo Prevete in formulating some of the R code above during review###################
########################################################################################################################