# Starting with ordinations ----

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

# Saving data in R

data(iris)
save(iris, file = "iris-df.rda")

iris <- "Iris is a beautiful flower."

load(file = "iris-df.rda")

saveRDS(iris, file = "iris-df.rds")

iris <- "Iris is a beautiful flower."

iris_the_table <- readRDS(file = "iris-df.rds")

## PCA ----
lalashan <- readRDS(file = "lalashan.rds")

names(lalashan)
head(lalashan$obs)

summary(lalashan$cross)

dim(lalashan$cross)

count_spp <- function(x) length(unique(x))

spp_freq <- aggregate(plot_id ~ taxon_code,
                      FUN = count_spp,
                      data = lalashan$obs)
head(spp_freq)

hist(spp_freq$plot_id)

## Displaying composition as dimensions
spp_dim <- lalashan$obs

# Sort the frequency of species
spp_freq <- spp_freq[order(spp_freq$plot_id, decreasing = TRUE),]
head(spp_freq)

subset(lalashan$spp, tolower(taxon_code) == "prunphae")

spp_dim$taxon_code <- factor(spp_dim$taxon_code,
                             levels = spp_freq$taxon_code)

# Graph species as dimensions
library(ggplot2)
library(viridis)

ggplot(data = spp_dim,
       aes(x = taxon_code, y = abund,
           group = plot_id, color = plot_id)) +
  geom_line() +
  geom_point() +
  scale_color_viridis(discrete = TRUE) +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 0.5, hjust = 1))

### All in one again! ----
library(ggplot2)
library(viridis)
lalashan <- readRDS(file = "lalashan.rds")
count_spp <- function(x) length(unique(x))
spp_freq <- aggregate(plot_id ~ taxon_code,
                      FUN = count_spp,
                      data = lalashan$obs)
spp_dim <- lalashan$obs
spp_freq <- spp_freq[order(spp_freq$plot_id, decreasing = TRUE),]
spp_dim$taxon_code <- factor(spp_dim$taxon_code,
                             levels = spp_freq$taxon_code)
ggplot(data = spp_dim,
       aes(x = taxon_code, y = abund,
           group = plot_id, color = plot_id)) +
  geom_line() +
  geom_point() +
  scale_color_viridis(discrete = TRUE) +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 0.5, hjust = 1))

## Principal Component Analysis ----
pca_ord <- prcomp(lalashan$cross)
biplot(pca_ord)

boxplot(abund ~ taxon_code, data = spp_dim)

abundance <- aggregate(abund ~ taxon_code,
                       data = spp_dim, FUN=mean)
abundance <- abundance[order(abundance$abund,
                             decreasing = TRUE),]
head(abundance)

# Run PCA again
pca_ord <- prcomp(lalashan$cross, scale. = TRUE)
biplot(pca_ord)

rare_spp <- spp_freq$taxon_code[spp_freq$plot_id < 5]
rare_spp

lalashan$cross2 <- lalashan$cross[,
                  !colnames(lalashan$cross) %in%
                    rare_spp]

pca_ord <- prcomp(lalashan$cross2, scale. = TRUE)
biplot(pca_ord)

# Show single species
lalashan$cross[,c("melisqua", "tricdubi")]
lalashan$cross[,c("euryglab", "querlong")]

# Exploring results
str(pca_ord)

plot(pca_ord)

biplot(pca_ord, choices = c(1,3))
biplot(pca_ord)

summary(pca_ord)