---
title: "Core Code for Ecological Quadrat Count Data" 
author: Chris 
date: "`r format(Sys.time(), '%d %B, %Y')`"\
output: 
  html_document: 
  theme: readable 
  highlight: zenburn 
    toc: true 
    toc_float: true 
    number_sections: true 
    code_folding: show 
    code_download: true 
editor: 
markdown: 
wrap: 72 
---

```{r setup, include=FALSE}
library(knitr)
install.packages("tidyverse")
library("tidyverse")
#library(kableExtra)
install.packages("vegan")
library(vegan)
install.packages("indicspecies")
library("indicspecies")
library(RColorBrewer)
library(data.table)
library('ggplot2')
knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = FALSE,         # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  fig.width = 6,       # Set plot width in inches
  fig.height = 4,      # Set plot height in inches
  fig.align = "center" # Align plots to the center
)
```

# General code for the basic analyses we will perform on the quadrat count data taken from yellow. When you use it, copy it to a new file so that you know what the original is and what your version is - this will help you to keep track of changes/ names/ etc.

#Setup

```{r}
# Load required libraries install.packages("vegan") # If not installed previously library(vegan)

# Read the data from CSV

intertidal_data <- read.csv("intertidal_data.csv", stringsAsFactors = FALSE)
```

# Species Richness

```{r}
species_richness <- length(unique(intertidal_data$Species))
```

# Species Diversity (Shannon's Diversity Index)

```{r}
diversity_index <- diversity(intertidal_data$Percent_Cover, index = "shannon")
```

# Species Abundance

```{r}
species_abundance <- aggregate(Quadrat_Count ~ Species, intertidal_data, sum)
```

# Community Composition visualization via PCA Plot

```{r}
# This is a plot that shows 3 dimensions of data in 2 dimensions
# Create a species abundance matrix

species_matrix <- t(table(intertidal_data$Species, intertidal_data$Quadrat_Count))

# Run PCA

pca_result <- rda(species_matrix)

# Create a PCA plot in base R - do this first and then do it in ggplot

plot(pca_result, type = "points", display = "species", cex = 1.5)

# Load required libraries

# install.packages("ggplot2") # If not installed previously library(ggplot2)

# Extract PCA scores for each species

pca_scores <- scores(pca_result, display = "species")

# Convert PCA scores to a data frame

pca_data <- as.data.frame(pca_scores$species)

# Add Species names as a column

pca_data$Species <- rownames(pca_data)

# Create a PCA plot using ggplot2

ggplot(pca_data, aes(x = PC1, y = PC2, label = Species)) + geom_point() + geom_text(size = 3, hjust = 0, vjust = 0) + xlab("Principal Component 1") + ylab("Principal Component 2") + ggtitle("PCA Plot of Intertidal Species") + theme_minimal()
```

#Species Accumulation Curve

```{r}
# Load required library install.packages("vegan") # If not installed previously library(vegan)

# Assuming you have a data frame 'species_data' with species names in column 'Species'. We may have to do this with genus or even family

# and quadrat IDs in column 'Quadrat_ID'

# Calculate species accumulation curve

accum_curve <- specaccum(table(species_data$Species, species_data$Quadrat_ID))

# Plot the species accumulation curve

plot(accum_curve, xlab = "Number of Quadrats Sampled", ylab = "Number of Species")
```

# Species - Environment Relationships

```{r}
# Load required library

# install.packages("vegan") # If not installed previously library(vegan)

# Assuming you have environmental variables in columns 'Temperature', 'Salinity', etc. - We Only have substrate type but we can use this later if we find environmental data from NOAA or the likes

# Create a community data matrix

community_matrix <- table(species_data$Species, species_data$Quadrat_ID)

# Create an environmental data matrix

environment_matrix <- as.matrix(species_data[, c("Temperature", "Salinity")])

# Run Canonical Correspondence Analysis (CCA)

cca_result <- cca(community_matrix ~ ., data = species_data[, c("Species", "Quadrat_ID")], env = environment_matrix)

# Summary of the CCA results

summary(cca_result)
```

# Cluster Analysis

```{r}
# Assuming you have a data frame 'species_data' with species abundances in quadrats/ This is used to compare between sites/ sections for us

# Calculate Bray-Curtis dissimilarity matrix

bray_curtis_matrix <- vegdist(species_data[, -c("Species", "Quadrat_ID")], method = "bray")

# Perform hierarchical clustering

hc <- hclust(bray_curtis_matrix)

# Plot the dendrogram

plot(hc, main = "Dendrogram of Intertidal Quadrats", xlab = "Quadrats", ylab = "Distance")
```

# Indicator Species Analysis

```{r}
# Assuming you have a data frame 'species_data' with species abundances and environmental data

# Load required library

# install.packages("indicspecies") # If not installed previously library(indicspecies)

# Run Indicator Species Analysis

indicator_species <- multipatt(species_data[, -c("Species", "Quadrat_ID")], species_data$Species)

# View the results

summary(indicator_species)
```

# Non-metric Multidimensional Scaling (NMDS)

```{r}
# Load required library

# install.packages("vegan") # If not installed previously library(vegan)

# Assuming you have a data frame 'species_data' with species abundances in quadrats

# Calculate Bray-Curtis dissimilarity matrix

bray_curtis_matrix <- vegdist(species_data[, -c("Species", "Quadrat_ID")], method = "bray")

# Perform NMDS

nmds_result <- metaMDS(bray_curtis_matrix, distance = "bray")

# Plot the NMDS ordination

plot(nmds_result, type = "text")
```

# Permutational Analysis of Variance (PERMANOVA)

```{r}

# Load required library

# install.packages("vegan") # If not installed previously library(vegan)

# Assuming you have a data frame 'species_data' with species abundances in quadrats

# Calculate Bray-Curtis dissimilarity matrix

bray_curtis_matrix <- vegdist(species_data[, -c("Species", "Quadrat_ID")], method = "bray")

# Perform PERMANOVA

permanova_result <- adonis(formula = bray_curtis_matrix ~ Grouping_Variable, data = species_data)

# View the results

summary(permanova_result)
```

# Correlation Analysis

```{r}
# Assuming you have a data frame 'species_data' with species abundances and environmental data

# Calculate the correlation matrix

correlation_matrix <- cor(species_data[, -c("Species", "Quadrat_ID")])

# View the correlation matrix

print(correlation_matrix)
```

# Species Abundance Distribution (SAD) Analysis

```{r}
# Assuming you have a data frame 'species_data' with species abundances in quadrats

# Calculate species abundance distribution

sad <- table(species_data$Species)

# Plot the SAD

plot(sad, xlab = "Species Abundance", ylab = "Frequency", main = "Species Abundance Distribution")
```

# Heatmap

```{r}
#install.packages("reshape2")
# Load required libraries
library(ggplot2)

# Prepare data for heatmap
heatmap_data <- intertidal_data[, -1]  # Remove the "Site" column
rownames(heatmap_data) <- intertidal_data$Site

# Create the heatmap
heatmap_plot <- ggplot(heatmap_data, aes(x = Species, y = Site, fill = Species)) +
  geom_tile() +
  theme_minimal() +
  scale_fill_manual(values = c("white", "blue"), labels = c("Absent", "Present")) +
  labs(title = "Presence/Absence Heatmap",
       x = "Species",
       y = "Site")

print(heatmap_plot)

```

# Stacked Bar Plot

```{r}
# Load required libraries
library(ggplot2)
library(reshape2)

# Melt the data for stacked barplot
melted_data <- melt(intertidal_data, id.vars = "Site")

# Create the stacked barplot
stacked_barplot <- ggplot(melted_data, aes(x = Site, fill = factor(value))) +
  geom_bar() +
  theme_minimal() +
  scale_fill_manual(values = c("blue", "white"), labels = c("Present", "Absent")) +
  labs(title = "Presence/Absence Stacked Barplot",
       x = "Site",
       y = "Count",
       fill = "Presence")

print(stacked_barplot)

```

# Stacked Area Plot

```{r}
# Load required libraries
library(ggplot2)
library(reshape2)

# Melt the data for stacked area plot
melted_data <- melt(intertidal_data, id.vars = "Site")

# Create the stacked area plot
stacked_area_plot <- ggplot(melted_data, aes(x = Site, y = value, fill = variable)) +
  geom_area() +
  theme_minimal() +
  labs(title = "Intertidal Percent Cover",
       x = "Site",
       y = "Percent Cover",
       fill = "Species")

print(stacked_area_plot)

```