# Load packages
library(vegan)
library(vegtable)

# Load data
lalashan <- readRDS("lalashan.rds")

# Calculate frequencies
spp_freq <- aggregate(plot_id ~ taxon_code,
                      FUN = function(x) length(unique(x)),
                      data = lalashan$obs)
spp_freq <- spp_freq[order(spp_freq$plot_id, decreasing = TRUE), ]
head(spp_freq)
summary(as.factor(spp_freq$plot_id))

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

# Transform abundance values
hist(lalashan$obs$abund)
hist(log(lalashan$obs$abund + 1))

lalashan$obs$abund_log <- log(lalashan$obs$abund + 1)

# Delete rare species and do cross table
cross_table <- crosstable(abund_log ~ taxon_code + plot_id,
                          FUN = mean,
                          data = subset(lalashan$obs, taxon_code %in%
                                          sel_spp),
                          na_to_zero = TRUE,
                          as_matrix = TRUE)

# Canonical Correspondence Analysis
head(lalashan$env)
lalashan$env$zone <- substr(lalashan$env$plot_id,
                            start = 1,
                            stop = 2)
summary(as.factor(lalashan$env$zone))

lalashan$env$wind <- substr(lalashan$env$plot_id,
                            start = 3,
                            stop = 3)
summary(as.factor(lalashan$env$wind))

# Decomposing aspect
lalashan$env$northness <- cos(lalashan$env$aspect*(pi/180))
lalashan$env$eastness <- sin(lalashan$env$aspect*(pi/180))

with(lalashan$env,
     plot(aspect[order(aspect)], northness[order(aspect)],
          type = "b"))
with(lalashan$env,
     plot(aspect[order(aspect)], eastness[order(aspect)],
          type = "b"))

# Canonical Correspondence Analysis
discard_vars <- c("plot_id","elevation_zone", "lat", "lon",
                  "aspect", "plot_size")
env_table <- lalashan$env[!names(lalashan$env) %in% discard_vars]

cca_ord <- cca(cross_table ~ ., data = env_table)
plot(cca_ord)
plot(cca_ord, display = c("wa", "cn"))

# Stepwise selection ----
sample(letters, size = 3)
sample(letters, size = 3)

set.seed(123)
sample(letters, size = 3)

# Two models
mod0 <- cca(cross_table ~ 1, data = env_table)
mod1 <- cca(cross_table ~ ., data = env_table)

mod <- ordistep(mod0, scope = formula(mod1), trace = FALSE,
                permutations = how(nperm = 499))
mod
mod$anova

mod <- ordistep(mod0, scope = formula(mod1), trace = TRUE,
                permutations = how(nperm = 499),
                direction = "forward")
mod$anova

# Final Selection
set.seed(123)
mod0 <- cca(cross_table ~ 1, data = env_table)
mod1 <- cca(cross_table ~ ., data = env_table)

mod <- ordistep(mod0, scope = formula(mod1), trace = FALSE,
                permutations = how(nperm = 499))
mod
mod$anova

# Final Ordination
cca_ord <- cca(cross_table ~ zone + ph_kcl + fe + total_n,
               data = env_table)
plot(cca_ord)

# ANOVA using permutation
anova(cca_ord, by = "terms")

plot(cca_ord, display = c("wa", "cn"))
ordisurf(cca_ord ~ mean_coldest,
         data = lalashan$env, col = "red",
         add = TRUE)
ordisurf(cca_ord ~ frosty_days,
         data = lalashan$env, col = "darkgreen",
         add = TRUE)

plot(cca_ord, display = c("wa"))
fitted_env <- envfit(cca_ord ~ .,
                     data = env_table)
plot(fitted_env)

# correlation matrix of variables
cor(env_table[!names(env_table) %in%
                c("zone", "wind")])

pca_env <- prcomp(env_table[!names(env_table) %in%
                           c("zone", "wind")],
               scale. = TRUE)
biplot(pca_env)

## Distance based RDA 
mod0 <- dbrda(cross_table ~ 1, data = env_table,
              distance = "bray")
mod1 <- dbrda(cross_table ~ ., data = env_table,
              distance = "bray")

mod <- ordistep(mod0, scope = formula(mod1),
                trace = FALSE,
                permutations = how(nperm = 499))
mod
mod$anova

dbrda_ord <- dbrda(cross_table ~ zone +
                     cn_ratio +
                     heatload + s + cu,
                   data = env_table)
plot (dbrda_ord)
anova(dbrda_ord, by = "terms")

# Procrustes ----
ord1 <- metaMDS(cross_table, distance = "bray")
plot(ord1, display = "sites", type = "t")

ord2 <- metaMDS(
  env_table[!names(env_table) %in%
              c("zone", "wind")],
                distance = "euclidean")
plot(ord2, display = "sites", type = "t")

comp_ords <- protest(ord2, ord1)
plot(comp_ords)
plot(comp_ords, kind = 2)
barplot(residuals(comp_ords))