# Downloading required files ----
download.file(
  url = "https://kamapu.gitlab.io/multivar/data/course-data.zip",
  destfile = "course-data.zip", method = "curl")
unzip("course-data.zip", overwrite = TRUE)
unlink("course-data.zip")

# Loading required packages
library(vegan)
library(vegtable)

# Load data
wetlands <- readRDS("wetlands.rds")
names(wetlands)

# Redundancy Analysis
head(wetlands$env)
rda_ord <- rda(wetlands$cross ~ groundwater + tmean +
                 tseason + prsum + prseason,
               data = wetlands$env)
rda_ord

plot(rda_ord)

# Canonical Correspondece Analysis
cca_ord <- cca(wetlands$cross ~ groundwater + tmean +
                 tseason + prsum + prseason,
               data = wetlands$env)
cca_ord

plot(cca_ord)

# Trying our best assessment ----
# - Taxonomic level?
# - Abundance?
head(wetlands$obs)
hist(wetlands$obs$cover_perc)
summary(wetlands$obs$cover_class)

A <- cut(wetlands$obs$cover_perc,
         breaks = c(0, 25, 50, 70, 100))
summary(A)

A <- cut(wetlands$obs$cover_perc,
         breaks = c(0, 25, 50, 70, 100),
         labels = FALSE)
summary(as.factor(A))

wetlands$obs$class2 <- cut(wetlands$obs$cover_perc,
                           breaks = c(0, 25, 50, 70, 100),
                           labels = FALSE)
hist(wetlands$obs$class2)

# Skip rare species
spp_freq <- aggregate(plot_id ~ taxon_code,
              FUN = function(x) length(unique(x)),
              data = wetlands$obs)
head(spp_freq)

spp_freq <- spp_freq[order(spp_freq$plot_id, decreasing = TRUE), ]
head(spp_freq)
tail(spp_freq)

sel_spp <- with(spp_freq, taxon_code[plot_id > 4])

cross_table <- crosstable(class2 ~ taxon_code + plot_id,
                FUN = max,
                data = subset(wetlands$obs, taxon_code %in% sel_spp),
                na_to_zero = TRUE,
                as_matrix = TRUE)

# which ordination? ----
dca_ord <- decorana(cross_table)
plot(dca_ord)
abline(v = -2, lty = "dashed", col = "red")
abline(v = 2, lty = "dashed", col = "red")

# Correspondence Anlysis?
ca_ord <- cca(cross_table)
plot(ca_ord)

fitted_model <- envfit(dca_ord ~ groundwater + tmean +
                         tseason + prsum + prseason,
                       data = wetlands$env)

plot(dca_ord)
plot(fitted_model)

# Canonical Correspondence Analysis
cca_ord <- cca(cross_table ~ groundwater + tmean +
                 tseason + prsum + prseason,
               data = wetlands$env)
plot(cca_ord)

# Correlations
cor(wetlands$env[c("groundwater", "tmean", "tseason",
                   "prsum", "prseason")])
pairs(wetlands$env[c("groundwater", "tmean", "tseason",
                   "prsum", "prseason")])

library(GGally)
ggpairs(wetlands$env[c("groundwater", "tmean", "tseason",
                     "prsum", "prseason")],
        lower = list(continuous = "smooth"))

plot(cca_ord)
cca_ord
str(cca_ord)

barplot(cca_ord$CCA$eig)
abline(h=mean(cca_ord$CCA$eig), lty = 2, col = "red")

plot(cca_ord, choices = c(1,3))

par(mfrow = c(3, 1))
plot(cca_ord)
plot(cca_ord, choices = c(1,3))
plot(cca_ord, choices = c(2, 3))

# Explained variance in percentage
var_table <- with(cca_ord$CCA, {
  expl_var = eig/sum(eig)*100
  cum_var = cumsum(expl_var)
  data.frame(expl_var, cum_var)
})

var_table

par(mfrow = c(1,1))
plot(cca_ord, display = c("wa", "bp"))
text(-6, 5, labels = paste0("Explained Variance = ",
    round(var_table["CCA2", "cum_var"], digits = 1), "%"), pos =4)

# Selecting explanatory variables
cca_ord <- cca(cross_table ~ groundwater + tmean +
                 prsum + tseason:prseason,
               data = wetlands$env)
plot(cca_ord)

cca_ord <- cca(cross_table ~ groundwater + tmean +
                 prsum + tseason,
               data = wetlands$env)
plot(cca_ord)
text(-4, 5, labels = paste0("Explained Variance = ",
            round(var_table["CCA2", "cum_var"], digits = 1), "%"),
     pos = 4)

png(filename = "cca-ordination.png", width = 1500, height = 1500,
    res = 200)

par(mar = c(5, 5, 1, 1))
plot(cca_ord)
text(-4, 5, labels = paste0("Explained Variance = ",
                            round(var_table["CCA2", "cum_var"],
                                  digits = 1), "%"),
     pos =4)

dev.off()

# db-RDA ----
dbrda_ord <- dbrda(cross_table ~ groundwater + tmean + prsum +
                     tseason + prseason,
                   data = wetlands$env)


attach(wetlands$env)
tmean
detach(wetlands$env)

with(wetlands$env, {
  tmean
})


dbrda_ord <- dbrda(cross_table ~ groundwater + tmean + prsum +
                     tseason + prseason,
                   data = wetlands$env, distance = "bray")
plot(dbrda_ord)

dbrda_ord2 <- dbrda(cross_table ~ groundwater + tmean + prsum +
                     tseason + prseason,
                   data = wetlands$env, distance = "jaccard")
plot(dbrda_ord2)

dbrda_ord <- dbrda(cross_table ~ groundwater + tmean + prsum +
                     tseason + prseason,
                   data = wetlands$env, distance = "man")
plot(dbrda_ord)

dbrda_ord

# Final Ordination
dbrda_ord <- dbrda(cross_table ~ groundwater + tmean + prsum +
                     tseason,
                   data = wetlands$env, distance = "bray")
plot(dbrda_ord)