LONG-TERM CORAL DATA FROM CPCE PHOTOQUAD ANNOTATIONS ############################
# load packages
library(readxl)
library(dplyr)
library(tidyverse)
library(vegan)
library(ggordiplots)
library(factoextra)
# load raw datasets
coral_2006 <- read_excel("Raw_Data/qry_benthic_data_all_2006.xlsx")
coral_2007 <- read_excel("Raw_Data/qry_benthic_data_all_2007.xlsx")
coral_2008 <- read_excel("Raw_Data/qry_benthic_data_all_2008.xlsx")
coral_2009 <- read_excel("Raw_Data/qry_benthic_data_all_2009.xlsx")
coral_2010 <- read_excel("Raw_Data/qry_benthic_data_all_2010.xlsx")
coral_2017 <- read_excel("Raw_Data/qry_benthic_data_all_2017.xlsx")
coral_2018 <- read_excel("Raw_Data/qry_benthic_data_all_2018.xlsx")
coral_2019 <- read_excel("Raw_Data/qry_benthic_data_all_2019.xlsx")# Function to process datasets
process_df <- function(df) {
# List of non-coral labels
non_coral <- c("Turf algae", "Coralline algae", "Unknown Coral", "Rubb", "Non-Coral",
"Sand","Blue green","Caulerpa serrulata", "Dictyota sp", "Unknown",
"Macroalgae", "Bare Rock","Echinoidea", "Asparagopsis sp","Shadow", "Acanthophora spicifera",
"Unindentified green", "Caulerpa racemosa", "Zoanthidae",
"Holothuriidae","Dictyosphaeria versluysii", "Padina sp", "Acanthaster planci" ,
"Holothuria sp", "Unindentified red", "Porifera","Wand","Tape","Black Hole","Gelatinous red",
"Sargassum sp","Bryozoa","Other sessile inverts","Liagora sp","Jointed calcareous red",
"Other Mobile Inverts","Galaxaura sp")
# Replace non-coral values with "Other"
df$Taxon_Name <- ifelse(df$Taxon_Name %in% non_coral, "Other", df$Taxon_Name)
# Compute % cover for each coral species
counts <- df %>%
group_by(Year, Loc_Name, Frame, Taxon_Name) %>%
summarise(count = n(), .groups = "drop") %>%
ungroup() %>%
group_by(Loc_Name, Frame, Year) %>%
mutate(total_count = sum(count)) %>%
mutate(cover_percent = ifelse(total_count == 0, 0, (count / total_count) * 100)) %>%
select(-count, -total_count)
# Compute overall % coral cover
counts <- counts %>%
group_by(Loc_Name, Frame) %>%
mutate(coral_cover = sum(cover_percent[!Taxon_Name == "Other"]))
# Pivoting so that each spp has its own column
df_final <- pivot_wider(counts, names_from = Taxon_Name, values_from = cover_percent,
values_fill = list(cover_percent = 0))
return(df_final)
}# Process each dataframe by applying the function above, and store them all in a list
cover_dataframes <- lapply(list(coral_2006, coral_2007, coral_2008, coral_2009,
coral_2010, coral_2017, coral_2018, coral_2019), process_df)
# Combine dataframes
KAHO_CoralCover <- bind_rows(cover_dataframes) %>% # combine each year dataframe into one
# make these variables a factor type variable
mutate_at(c("Year","Loc_Name","Frame"), factor) %>%
# replace NAs with 0
mutate(across(where(is.numeric), ~ifelse(is.na(.), 0, .))) %>%
# Specify the order wanted for the Year variable
mutate(Year = factor(as.character(Year),
levels = c("2006", "2007", "2008", "2009", "2010", "2017", "2018", "2019"))) %>%
ungroup()
# Summarize the data by computing mean cover per transect ("Loc_Name")
CoralCover_summ <- KAHO_CoralCover %>%
group_by(Year, Loc_Name) %>%
summarize(across(where(is.numeric), mean, na.rm = TRUE)) # take mean by "Loc_Name" for each "Year"Visualize Coral Cover Across Years
ggplot(CoralCover_summ, aes(x = Year, y = coral_cover, fill = Year)) +
geom_boxplot() + theme_classic()
ggplot(CoralCover_summ, aes(x = Year, y = `Porites lobata`, fill = Year)) +
geom_boxplot() + theme_classic()
ggplot(CoralCover_summ, aes(x = Year, y = `Pocillopora meandrina`, fill = Year)) +
geom_boxplot() + theme_classic()
ggplot(CoralCover_summ, aes(x = Year, y = `Porites compressa`, fill = Year)) +
geom_boxplot() + theme_classic()
ggplot(CoralCover_summ, aes(x = Year, y = `Porites evermanni`, fill = Year)) +
geom_boxplot() + theme_classic()
ggplot(CoralCover_summ, aes(x = Year, y = `Montipora capitata`, fill = Year)) +
geom_boxplot() + theme_classic()
Mean % Coral Cover by Year with Standard Error bars and Trend Line
CoralCover_summ2 <- KAHO_CoralCover %>%
group_by(Year) %>%
summarize(mean_coral_cover = mean(coral_cover, na.rm = TRUE),
se_coral_cover = sd(coral_cover, na.rm = TRUE) / sqrt(n()))
ggplot(CoralCover_summ2, aes(x = Year, y = mean_coral_cover, group = 1)) +
geom_line() + geom_point(size = 1) +
geom_errorbar(aes(ymin = mean_coral_cover - se_coral_cover,
ymax = mean_coral_cover + se_coral_cover), width = 0.2) +
theme_bw()
Non-Metric Multidimensional Scaling - Coral Community Across Years
KAHO_ord <- KAHO_CoralCover %>% dplyr::select(-c("Frame", "Other","coral_cover", "Loc_Name"))
# clean up to common coral because there's an issue with distance calculations
KAHO_ord <- KAHO_ord %>%
dplyr::select(-c("Octocoral", "Porites massive sp","Acanthophora sp", "Lobophyllia sp",
"Sarcothelia edmondsoni", "Porites sp","Leptoseris incrustans",
"Leptastrea bewickensis","Leptastrea purpurea","Cyphastrea ocellina","Pavona duerdeni"))
# remove sites where there was % 0 coral cover
KAHO_ord <- KAHO_ord[ rowSums(KAHO_ord[,-1]) > 0, ]
# remove outliers (spotted on first NMDS run)
KAHO_ord <- KAHO_ord %>%
slice(-2580)
# remove outliers (spotted on second NMDS run)
KAHO_ord <- KAHO_ord %>%
slice(-2942)
# remove outliers (spotted on third NMDS run)
KAHO_ord <- KAHO_ord %>%
slice(-c(3501,3037))
# bray-curtid dissimilarity matrix on coral community data
KAHO_dist <- as.matrix(vegdist(KAHO_ord[,-1], method = "bray"))
# run NMDS
KAHO_nmds <- metaMDS(KAHO_dist, k=2)## Run 0 stress 0.1553449
## Run 1 stress 0.1598859
## Run 2 stress 0.1963284
## Run 3 stress 0.162791
## Run 4 stress 0.1658176
## Run 5 stress 0.1637609
## Run 6 stress 0.1641969
## Run 7 stress 0.1580466
## Run 8 stress 0.158826
## Run 9 stress 0.1665643
## Run 10 stress 0.1705642
## Run 11 stress 0.1635499
## Run 12 stress 0.1831204
## Run 13 stress 0.1714078
## Run 14 stress 0.1651366
## Run 15 stress 0.1616716
## Run 16 stress 0.1622408
## Run 17 stress 0.1813402
## Run 18 stress 0.173884
## Run 19 stress 0.1696516
## Run 20 stress 0.1673911
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 9: no. of iterations >= maxit
## 11: scale factor of the gradient < sfgrmin# Shepard diagram = ordination distances and linear fit line against original dissimilarities
stressplot(KAHO_nmds)
# fit species vectors to the ordination
en = envfit(KAHO_nmds, KAHO_ord[,-1], permutations = 999, na.rm = TRUE)
# extract NMDS scores data
scrs <- as.data.frame(scores(KAHO_nmds, display = "sites"))
spp.scrs <- as.data.frame(scores(en, display = "vectors"))
spp.scrs <- cbind(spp.scrs, Species = rownames(spp.scrs))
ggplot(scrs) +
# add points from the first two NMDS axes and colour them by Year
geom_point(mapping = aes(x = NMDS1, y = NMDS2, colour = KAHO_ord$Year), size =1) +
# add ellipse with 95% CI over each year's data
stat_ellipse(mapping = aes(x = NMDS1, y = NMDS2, colour = KAHO_ord$Year),
level = 0.95, type = "norm") +
coord_fixed() + # set fixed aspect ratio
# add species vectors
geom_segment(data = spp.scrs,
aes(x = 0, xend = 0.4 * NMDS1, y = 0, yend = 0.4 * NMDS2),
arrow = arrow(length = unit(0.25, "cm")), colour = "black") +
# add species labels
geom_text(data = spp.scrs, aes(x = 0.43 * NMDS1, y = 0.43 * NMDS2, label = Species),
size = 4.6) + theme_classic() +
theme(legend.position = "bottom", # specify legend position
legend.title = element_text(size = 12, face = "bold"), # adjust legend title size and style
legend.text = element_text(size = 10), # adjust legend text size
legend.key.size = unit(4, "mm")) + # adjust legend symbol size
scale_colour_discrete(name = "Year") # adjust legend title
# PCA
spp_pca <- prcomp(KAHO_ord[,-1], scale = FALSE)
factoextra::fviz_pca_var(spp_pca,col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel = TRUE)