# DCA ----
# Load the data

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")

# Lalashan
lalashan <- readRDS("lalashan.rds")

# Load the package
library(vegan)

dca <- decorana(lalashan$cross)
dca

plot(dca)

head(lalashan$obs)
boxplot(abund ~ taxon_code, data = lalashan$obs)

# standardizing
?scale
cross_stand <- scale(lalashan$cross)

lalashan$cross[1:10,1:10]
cross_stand[1:10,1:10]

dca <- decorana(cross_stand)

# Take a look on the abundance
summary(lalashan$cross[,1:10])
summary(cross_stand[,1:10])

apply(lalashan$cross, 2, FUN = mean)
apply(cross_stand, 2, FUN = function(x) round(mean(x)))
apply(cross_stand, 2, FUN = function(x) round(sd(x)))

cross_stand <- scale(lalashan$cross, center = FALSE)

dca <- decorana(cross_stand)
plot(dca)

# Calculate frequencies
spp_freq <- aggregate(plot_id ~ taxon_code,
                      data = lalashan$obs,
                      FUN = function(x) length(unique(x)))
head(spp_freq)

sel_spp <- with(spp_freq, {
  taxon_code[plot_id > 3]
})

sel_spp <- spp_freq$taxon_code[spp_freq$plot_id > 3]

cross2 <- subset(lalashan$obs, taxon_code %in% sel_spp)
head(cross2)

library(vegtable)
cross2 <- crosstable(abund ~ taxon_code + plot_id,
                     data = cross2,
                     FUN = sum,
                     na_to_zero = TRUE,
                     as_matrix = TRUE)

cross2 <- scale(cross2, center = FALSE)

dca <- decorana(cross2)
dca

plot(dca)

# Some exploration
lalashan$spp["michcomp",]

lalashan$cross[,"michcomp"]
lalashan$cross[,c("michcomp","ilexfico")]

# Compressed code
spp_freq <- aggregate(plot_id ~ taxon_code,
                      data = lalashan$obs,
                      FUN = function(x) length(unique(x)))
sel_spp <- spp_freq$taxon_code[spp_freq$plot_id > 3]
cross2 <- subset(lalashan$obs, taxon_code %in% sel_spp)
library(vegtable)
cross2 <- crosstable(abund ~ taxon_code + plot_id,
                     data = cross2,
                     FUN = sum,
                     na_to_zero = TRUE,
                     as_matrix = TRUE)
cross2 <- scale(cross2, center = FALSE)
dca <- decorana(cross2)
plot(dca)

vignette(package="vegan")
vignette("intro-vegan")

# Use overlays (post-hoc correlations) ----
plot_ordfit <- envfit(dca ~ elevation + frosty_days,
                      data = lalashan$env)
plot(plot_ordfit)

ordisurf(dca ~ elevation, data = lalashan$env, add = TRUE)

str(dca)

# Customize plot
plot(dca, type = "n")
points(dca, display = "sites", pch = 21)
text(dca, display = "sites", pos = 4)

# Plot again
plot(dca)
abline(v=c(-1.5, 1.5), lty = 2, col = "red")
abline(v=c(-2, 2), lty = 2, col = "darkgreen")

plot(dca)
plot(plot_ordfit)

plot(lalashan$env$elevation, lalashan$env$frosty_days)
with(lalashan$env, plot(elevation, frosty_days))

plot(dca)
plot(plot_ordfit)

# Non-Metric Multidimensional Scaling ----

wetlands <- readRDS("wetlands.rds")
head(Wetlands$ob)

library(vegan)
?vegdist

d_euc <- vegdist(wetlands$cross, method = "euclidean")

as.matrix(d_euc)[1:10,1:10]
summary(d_euc)

d_bray <- vegdist(wetlands$cross, method = "bray")
as.matrix(d_bray)[1:10,1:10]
summary(d_bray)

nmds_ord <- metaMDS(comm = wetlands$cross,
                    distance = "bray",
                    k = 2)
plot(nmds_ord)