Exploratory Ecological Analysis of MERMAID Data – Rare Honduras

Summary

In this document, I explore the available fish and benthic monitoring data for Rare Honduras, all of which is currently uploaded and accessible through the MERMAID platform. The analysis includes a set of initial exploratory visualizations intended to accompany the Ecological Report and Actionable Recommendations for Managed Access & Reserves in Honduras currently being prepared for BAF. Additional visualizations can be produced based on feedback or specific needs identified by the team following this initial review.

---

Overview of Available Ecological Data

# Get MERMAID Honduras Projects (sites and year):
# Authenticate (use your MERMAID credentials or token)
mermaid_auth()

# Get all projects, including test projects
all_projects <- mermaid_get_projects(include_test_projects = TRUE)

# Filter for Rare Honduras projects
rare_projects_honduras <- all_projects %>%
  filter(grepl("Rare", tags, ignore.case = TRUE), countries == "Honduras") %>%
  select(project_id = id, project_name = name)


# Define methods to check
methods <- c("fishbelt", "benthiclit", "benthicpit", "benthicpqt", "bleaching", "habitatcomplexity")

# Initialize results list
results <- list()

# Loop through the first 10 projects
for (i in 1:10) {
  project_id <- rare_projects_honduras$project_id[i]
  project_name <- rare_projects_honduras$project_name[i]
  
  method <- NA
  
  for (m in methods) {
    obs <- tryCatch(
      mermaid_get_project_data(project = project_id, method = m, data = "observations"),
      error = function(e) NULL
    )
    
    if (is.data.frame(obs) && nrow(obs) > 0) {
      method <- m
      break
    }
  }
  
  results[[i]] <- data.frame(
    project_id = project_id,
    project_name = project_name,
    method = method,
    stringsAsFactors = FALSE
  )
}

# Combine into one data frame
project_methods <- bind_rows(results)

# ⚠️ Manually fix method for known Benthic Photo Quadrat projects
# This is necessary because MERMAID's API does not return data for benthicpqt!!!
project_methods$method[project_methods$project_name %in% c(
  "honduras_coralnet_2021-2023_habitat",
  "honduras_ff2.0_habitat"
)] <- "benthicpqt"


# Define the managed access areas in the same order as the projects
managed_access_areas <- c(
  "Puerto Cortes",
  "Santa Fe",
  "Guanaja; Iriona and Limon; Puerto Cortes; Roatan; Santa Fe; Trujillo; Utila",
  "Trujillo",
  "Trujillo; Santa Fe; Guanaja; Roatan; Iriona and Limon",
  "Belfate",
  "Guanaja",
  "Utila",
  "Iriona and Limon; Puerto Cortes; Roatan; Santa Fe; Utila",
  "Roatan"
)

# Add the column to the data frame
project_methods$managed_access_area <- managed_access_areas


# Clean and select relevant columns
project_methods_clean <- project_methods %>%
  select(
    `Project Name` = project_name,
    `Method` = method,
    `Managed Access Area` = managed_access_area
  )

# Generate the table
project_methods_clean %>%
  kbl(
    caption = "Projects currently available in MERMAID from Rare in Honduras",
    align = "l"
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Projects currently available in MERMAID from Rare in Honduras

Project Name	Method	Managed Access Area
Ecological Monitoring Bahía Puerto Cortés	fishbelt	Puerto Cortes
SantaFe_2022	fishbelt	Santa Fe
honduras_ff2.0_fish	fishbelt	Guanaja; Iriona and Limon; Puerto Cortes; Roatan; Santa Fe; Trujillo; Utila
Ecological Monitoring ZRP Bajo Calderon	fishbelt	Trujillo
honduras_coralnet_2021-2023_habitat	benthicpqt	Trujillo; Santa Fe; Guanaja; Roatan; Iriona and Limon
Balfate Fish Survey, May 2024	fishbelt	Belfate
Guanaja_2022	fishbelt	Guanaja
Ecological Monitoring Utila (2 Reserves)	fishbelt	Utila
honduras_ff2.0_habitat	benthicpqt	Iriona and Limon; Puerto Cortes; Roatan; Santa Fe; Utila
Ecological Monitoring Roatan (2 Reserves)	fishbelt	Roatan

### Fish Data

# Getting all fishbelt data from MERMAID

# # Basic template for extracting MERMAID info for method = "fishbelt"
# project_id <- rare_projects_honduras$project_id[2]
# 
# # Retrieve sample event data for the selected project
# project_sample_events <- mermaid_get_project_data(
#   project = project_id,
#   method = "fishbelt",      # Replace with the appropriate method if different
#   data = "sampleevents"     # Options include "sampleevents", "sampleunits", "observations"
# )
# 
# 
# # Retrieve sample units data for the selected project
# project_sample_units <- mermaid_get_project_data(
#   project = project_id,
#   method = method,      # Replace with the appropriate method if different
#   data = "sampleunits"     # Options include "sampleevents", "sampleunits", "observations"
# )
# 
# 
# # Retrieve observations data for the selected project
# project_sample_observations <- mermaid_get_project_data(
#   project = project_id,
#   method = method,      # Replace with the appropriate method if different
#   data = "observations"     # Options include "sampleevents", "sampleunits", "observations"
# )


# Project indexes you want to include
selected_indexes <- c(1, 2, 3, 4, 6, 7, 8, 10)
selected_projects <- rare_projects_honduras[selected_indexes, ]

# 1. Get sample events
sample_events_list <- map(selected_indexes, function(i) {
  pid <- rare_projects_honduras$project_id[i]
  pname <- rare_projects_honduras$project_name[i]
  
  tryCatch({
    mermaid_get_project_data(project = pid, method = "fishbelt", data = "sampleevents") %>%
      select(project = project, site, management, sample_date, biomass_kgha_avg) %>%
      mutate(project = pname)
  }, error = function(e) NULL)
})

# Combine all into one df
sample_events_df <- bind_rows(sample_events_list)



# 2. Get sample units
sample_units_list <- map(selected_indexes, function(i) {
  pid <- rare_projects_honduras$project_id[i]
  pname <- rare_projects_honduras$project_name[i]
  
  tryCatch({
    mermaid_get_project_data(project = pid, method = "fishbelt", data = "sampleunits") %>%
      select(project = project, site, management, sample_date, biomass_kgha, total_abundance) %>%
      mutate(project = pname)
  }, error = function(e) NULL)
})

# Combine all into one df
sample_units_df <- bind_rows(sample_units_list)


# 3. Get sample observations
sample_observations_list <- map(selected_indexes, function(i) {
  pid <- rare_projects_honduras$project_id[i]
  pname <- rare_projects_honduras$project_name[i]
  
  tryCatch({
    mermaid_get_project_data(project = pid, method = "fishbelt", data = "observations") %>%
      mutate(project = pname)
  }, error = function(e) NULL)
})

# Combine into a single dataframe
sample_observations_df <- bind_rows(sample_observations_list)




# # View the outputs
# View(sample_events_df)
# View(sample_units_df)
# View(sample_observations_df)


# Replace site names
sample_events_df <- sample_events_df %>%
  mutate(`site` = str_replace_all(str_to_lower(`site`), " ", "_"))
sample_units_df <- sample_units_df %>%
  mutate(`site` = str_replace_all(str_to_lower(`site`), " ", "_"))
sample_observations_df <- sample_observations_df %>%
  mutate(`site` = str_replace_all(str_to_lower(`site`), " ", "_"))

# Add MAA Names
sites <- read_excel(here("data", "processed", "sites", "honduras", "hon_ff2.0_fish_sites_processed.xlsx"))


# Perform a left join to bring in managed access area names
sample_events_df <- sample_events_df %>%
  left_join(sites %>% select(name, ma_name), by = c("site" = "name")) %>%
  rename(managed_access_areas = ma_name) 

sample_events_df <- sample_events_df %>%
  mutate(managed_access_areas = case_when(
    is.na(managed_access_areas) & project == "Balfate Fish Survey, May 2024" ~ "Belfate",
    is.na(managed_access_areas) & project == "Ecological Monitoring ZRP Bajo Calderon" ~ "Trujillo",
    is.na(managed_access_areas) & project == "Ecological Monitoring Roatan (2 Reserves)" ~ "Roatan",
    TRUE ~ managed_access_areas
  ))


sample_units_df <- sample_units_df %>%
  left_join(sites %>% select(name, ma_name), by = c("site" = "name")) %>%
  rename(managed_access_areas = ma_name) 

sample_units_df <- sample_units_df %>%
  mutate(managed_access_areas = case_when(
    is.na(managed_access_areas) & project == "Balfate Fish Survey, May 2024" ~ "Belfate",
    is.na(managed_access_areas) & project == "Ecological Monitoring ZRP Bajo Calderon" ~ "Trujillo",
    is.na(managed_access_areas) & project == "Ecological Monitoring Roatan (2 Reserves)" ~ "Roatan",
    TRUE ~ managed_access_areas
  ))


sample_observations_df <- sample_observations_df %>%
  left_join(sites %>% select(name, ma_name), by = c("site" = "name")) %>%
  rename(managed_access_areas = ma_name) 

sample_observations_df <- sample_observations_df %>%
  mutate(managed_access_areas = case_when(
    is.na(managed_access_areas) & project == "Balfate Fish Survey, May 2024" ~ "Belfate",
    is.na(managed_access_areas) & project == "Ecological Monitoring ZRP Bajo Calderon" ~ "Trujillo",
    is.na(managed_access_areas) & project == "Ecological Monitoring Roatan (2 Reserves)" ~ "Roatan",
    TRUE ~ managed_access_areas
  ))


###########################################
# Check for duplicates
# Check for duplicates based on site + sample_date + management
# 1. Normalize management field
sample_events_df <- sample_events_df %>%
  mutate(
    management = case_when(
      management %in% c("Managed Access", "Managed Access Area") ~ "Managed Access",
      TRUE ~ management
    )
  )

# 2. Identify duplicates based on site, sample_date, and cleaned management
fish_duplicated_events <- sample_events_df %>%
  group_by(site, sample_date, management) %>%
  filter(n() > 1) %>%
  arrange(site, sample_date, management)

# 3. View result
# nrow(fish_duplicated_events)         # Number of duplicate entries
# fish_duplicated_events               # Inspect duplicates


# Check for duplicates based on biomass + year from sample_date + maa + management:
# 1. Prepare the data by creating comparison columns
sample_events_df_with_flags <- sample_events_df %>%
  mutate(
    biomass_rounded = round(biomass_kgha_avg),
    year = year(sample_date)
  )

# 2. Identify duplicates based on rounded biomass, year, and MAA
biomass_duplicated_events <- sample_events_df_with_flags %>%
  group_by(biomass_rounded, year, managed_access_areas, management) %>%
  filter(n() > 1) %>%
  arrange(managed_access_areas, year, biomass_rounded)

# 3. View result
# nrow(biomass_duplicated_events)         # Number of duplicate entries
# biomass_duplicated_events               # Inspect duplicates



#Remove duplicates
# 1. Normalize management values
normalize_management <- function(df) {
  df %>%
    mutate(management = case_when(
      management %in% c("Managed Access", "Managed Access Area") ~ "Managed Access",
      TRUE ~ management
    ))
}

sample_events_df <- normalize_management(sample_events_df)
sample_units_df <- normalize_management(sample_units_df)
sample_observations_df <- normalize_management(sample_observations_df)

# 2. Create consistent site_key: site + sample_date + management
create_site_key <- function(df) {
  df %>%
    mutate(site_key = paste(site, sample_date, management, sep = "_"))
}

sample_events_df <- create_site_key(sample_events_df)
sample_units_df <- create_site_key(sample_units_df)
sample_observations_df <- create_site_key(sample_observations_df)

# 3. Identify site_keys that appear in more than one project in sample_events_df
site_key_projects <- sample_events_df %>%
  distinct(site_key, project) %>%
  count(site_key) %>%
  filter(n > 1) %>%
  pull(site_key)

# 4. Remove "honduras_ff2.0_fish" entries matching those site_keys
sample_events_df <- sample_events_df %>%
  filter(!(site_key %in% site_key_projects & project == "honduras_ff2.0_fish")) %>%
  select(-site_key)

sample_units_df <- sample_units_df %>%
  filter(!(site_key %in% site_key_projects & project == "honduras_ff2.0_fish")) %>%
  select(-site_key)

sample_observations_df <- sample_observations_df %>%
  filter(!(site_key %in% site_key_projects & project == "honduras_ff2.0_fish")) %>%
  select(-site_key)

# This code removes duplicate sampling events based on shared site name, sample date, and management type, specifically when the same event appears in multiple projects. Disclaimer: This method relies on exact matching of site name and date. If sites were recorded under slightly different names (e.g., "punto_10" vs. "10") or if dates were not entered consistently, some duplicates may not be detected and could still remain in the data. 

### Benthic Data

# Getting coralnet habitat data from MERMAID (not possible to directly get that from MERMAID, had to download and process the data)

# Step 1: Load data
coralnet <- read_excel(here("data", "raw", "honduras_coralnet_2021_2023_habitat_benthicpqt_25_06_16.xlsx"), 
                 sheet = "Benthic PQT Obs")

# Step 2: Define mapping function for benthic attributes
recode_categories <- function(attribute) {
  case_when(
    # HARD CORAL
    attribute %in% c(
  "Hard coral", "Agaricia", "Agaricia tenuifolia", "Siderastrea", "Siderastrea siderea",
  "Porites", "Porites astreoides", "Porites porites", "Porites furcata", "Montastraea",
  "Dendrogyra cylindrus", "Orbicella", "Orbicella annularis", "Orbicella faveolata",
  "Acropora cervicornis", "Acropora palmata", "Pseudodiploria strigosa",
  "Colpophyllia natans", "Diploria labyrinthiformis", "Meandrina meandrites",
  "Pocillopora ligulata", "Millepora", "Millepora alcicornis", "Millepora complanata"
) ~ "Hard coral",
    
    # SOFT CORAL
    attribute %in% c(
      "Soft coral", "Gorgonia", "Plexaura", "Pseudopterogorgia", "Pterogorgia", "Eunicea",
      "Pseudoplexaura"
    ) ~ "Soft coral",
    
    # TURF ALGAE
    attribute %in% c("Turf algae") ~ "Turf algae",
    
    # MACROALGAE
    attribute %in% c("Macroalgae", "Dictyota", "Lobophora", "Peyssonnelia", "Green algae") ~ "Macroalgae",
    
    # SAND
    attribute %in% c("Sand", "Silt") ~ "Sand",
    
    # RUBBLE
    attribute %in% c("Rubble") ~ "Rubble",
    
    # CRUSTOSE CORALLINE ALGAE
    attribute %in% c("Crustose coralline algae") ~ "Other",
    
    # BARE SUBSTRATE
    attribute %in% c("Bare substrate") ~ "Bare substrate",
    
    # OTHER INVERTEBRATES
    attribute %in% c(
      "Cliona","Erythropodium caribaeorum", "Palythoa", "Palythoa caribaeorum", "Sponge",
      "Other invertebrates"
    ) ~ "Other invertebrates",
    
    # SEAGRASS
    attribute %in% c("Thalassia testudinum", "Syringodium filiforme") ~ "Other",
    
    # UNKNOWN
    attribute %in% c("Unknown") ~ "Other",
    
    # DEFAULT
    TRUE ~ "Other"
  )
}

# Step 3: Apply mapping and calculate % cover by site + transect
# Step 3: Apply mapping and calculate % cover by site + transect + date
benthic_summary_coralnet <- coralnet %>%
  mutate(
    site = as.character(Site),
    transect = `Transect number`,
    sample_date = as.Date(paste(Year, Month, Day, sep = "-")),
    category = recode_categories(`Benthic attribute`)
  ) %>%
  filter(!is.na(site), !is.na(transect), !is.na(category)) %>%
  group_by(site, transect, sample_date, category) %>%
  summarise(point_count = sum(`Number of points`, na.rm = TRUE), .groups = "drop") %>%
  group_by(site, transect, sample_date) %>%
  mutate(
    total_points = sum(point_count),
    percent_cover = round(100 * point_count / total_points, 1)
  ) %>%
  ungroup() %>%
  select(site, transect, sample_date, category, percent_cover)

# View a specific site/transect to compare with MERMAID UI
# filter(benthic_summary_coralnet, site == 11, transect == 1)



# Add management and date
# Step 1: Create a lookup table for site/transect metadata
site_transect_info <- coralnet %>%
  mutate(
    site = as.character(Site),
    transect = `Transect number`,
    sample_date = as.Date(paste(Year, Month, Day, sep = "-"))
  ) %>%
  select(site, transect, management = `Management name`, sample_date) %>%
  distinct()

# Step 2: Join to benthic_summary
benthic_summary_coralnet <- benthic_summary_coralnet %>%
  left_join(site_transect_info, by = c("site", "transect", "sample_date")) %>%
  select(site, transect, sample_date, management, category, percent_cover)



# Add MAA Names
# Step 1: Read all site files
sites_coralnet_1 <- read_excel(here("data", "processed", "sites", "honduras", "hon_coralnet_arrecifes_habitat_sites_processed.xlsx"))
sites_coralnet_2 <- read_excel(here("data", "processed", "sites", "honduras", "hon_coralnet_calabash_habitat_sites_processed.xlsx"))
sites_coralnet_3 <- read_excel(here("data", "processed", "sites", "honduras", "hon_coralnet_guanaja_habitat_sites_processed.xlsx"))
sites_coralnet_4 <- read_excel(here("data", "processed", "sites", "honduras", "hon_coralnet_iriona_habitat_sites_processed.xlsx"))

# Step 2: Bind all into one lookup table
sites_all <- bind_rows(sites_coralnet_1, sites_coralnet_2, sites_coralnet_3, sites_coralnet_4)


# Remove duplicated benthic attributes per site, transect, date, and management
benthic_summary_coralnet <- benthic_summary_coralnet %>%
  group_by(site, transect, sample_date, management, category) %>%
  summarise(percent_cover = first(percent_cover), .groups = "drop")

# Step 3: Join to your summary table
benthic_summary_coralnet <- benthic_summary_coralnet %>%
  left_join(sites_all %>% select(name, site_name), by = c("site" = "name")) %>%
  rename(managed_access_areas = site_name)

#Remove tildes
benthic_summary_coralnet <- benthic_summary_coralnet %>%
  mutate(managed_access_areas = stri_trans_general(managed_access_areas, "Latin-ASCII"))


###########################################

# Getting FF 2.0 habitat data from MERMAID (not possible to directly get that from MERMAID, had to download and process the data)

# Step 1: Load data
ff2_habitat <- read_excel(here("data", "raw", "honduras_ff20_habitat_benthicpqt_25_06_16.xlsx"), 
                 sheet = "Benthic PQT Obs")

# Step 2: Define mapping function for benthic attributes
recode_categories <- function(attribute) {
  case_when(
    # HARD CORAL
    attribute %in% c(
  "Hard coral", "Agaricia", "Agaricia tenuifolia", "Siderastrea", "Siderastrea siderea",
  "Porites", "Porites astreoides", "Porites porites", "Porites furcata", "Montastraea",
  "Dendrogyra cylindrus", "Orbicella", "Orbicella annularis", "Orbicella faveolata",
  "Acropora cervicornis", "Acropora palmata", "Pseudodiploria strigosa",
  "Colpophyllia natans", "Diploria labyrinthiformis", "Meandrina meandrites",
  "Pocillopora ligulata", "Millepora", "Millepora alcicornis", "Millepora complanata"
) ~ "Hard coral",
    
    # SOFT CORAL
    attribute %in% c(
      "Soft coral", "Gorgonia", "Plexaura", "Pseudopterogorgia", "Pterogorgia", "Eunicea",
      "Pseudoplexaura"
    ) ~ "Soft coral",
    
    # TURF ALGAE
    attribute %in% c("Turf algae") ~ "Turf algae",
    
    # MACROALGAE
    attribute %in% c("Macroalgae", "Dictyota", "Lobophora", "Peyssonnelia", "Green algae") ~ "Macroalgae",
    
    # SAND
    attribute %in% c("Sand", "Silt") ~ "Sand",
    
    # RUBBLE
    attribute %in% c("Rubble") ~ "Rubble",
    
    # CRUSTOSE CORALLINE ALGAE
    attribute %in% c("Crustose coralline algae") ~ "Other",
    
    # BARE SUBSTRATE
    attribute %in% c("Bare substrate") ~ "Bare substrate",
    
    # OTHER INVERTEBRATES
    attribute %in% c(
      "Cliona","Erythropodium caribaeorum", "Palythoa", "Palythoa caribaeorum", "Sponge",
      "Other invertebrates"
    ) ~ "Other invertebrates",
    
    # SEAGRASS
    attribute %in% c("Thalassia testudinum", "Syringodium filiforme") ~ "Other",
    
    # UNKNOWN
    attribute %in% c("Unknown") ~ "Other",
    
    # DEFAULT
    TRUE ~ "Other"
  )
}

# Step 3: Apply mapping and calculate % cover by site + transect
benthic_summary_ff2 <- ff2_habitat %>%
  mutate(
    site = as.character(Site),
    transect = `Transect number`,
    category = recode_categories(`Benthic attribute`)
  ) %>%
  filter(!is.na(site), !is.na(transect), !is.na(category)) %>%
  group_by(site, transect, category) %>%
  summarise(point_count = sum(`Number of points`, na.rm = TRUE), .groups = "drop") %>%
  group_by(site, transect) %>%
  mutate(
    total_points = sum(point_count),
    percent_cover = round(100 * point_count / total_points, 1)
  ) %>%
  ungroup() %>%
  select(site, transect, category, percent_cover)

# View a specific site/transect to compare with MERMAID UI
# filter(benthic_summary_ff2, site == 11, transect == 1)



# Add management and date
# Step 1: Create a lookup table for site/transect metadata
site_transect_info <- ff2_habitat %>%
  mutate(
    site = as.character(Site),
    transect = `Transect number`,
    sample_date = as.Date(paste(Year, Month, Day, sep = "-"))
  ) %>%
  select(site, transect, management = `Management name`, sample_date) %>%
  distinct()

# Step 2: Join to benthic_summary
benthic_summary_ff2 <- benthic_summary_ff2 %>%
  left_join(site_transect_info, by = c("site", "transect")) %>%
  select(site, transect, sample_date, management, category, percent_cover)


# Add MAA Names
# Step 1: Read all site files
sites_ff2 <- read_excel(here("data", "processed", "sites", "honduras", "hon_ff2.0_habitat_sites_processed.xlsx"))

# Step 3: Join to your summary table
benthic_summary_ff2 <- benthic_summary_ff2 %>%
  left_join(sites_ff2 %>% select(name, site_name), by = c("site" = "name")) %>%
  rename(managed_access_areas = site_name)

#Remove tildes
benthic_summary_ff2 <- benthic_summary_ff2 %>%
  mutate(managed_access_areas = stri_trans_general(managed_access_areas, "Latin-ASCII"))


###########################################

# Combining all Benthic data (adding source column)
benthic_summary_coralnet <- benthic_summary_coralnet %>% mutate(source = "coralnet")
benthic_summary_ff2 <- benthic_summary_ff2 %>% mutate(source = "ff2")

benthic_summary_combined <- bind_rows(benthic_summary_coralnet, benthic_summary_ff2)


# write.csv(benthic_summary_coralnet, "benthic_summary_coralnet.csv", row.names = FALSE)
# write.csv(benthic_summary_ff2, "benthic_summary_ff2.csv", row.names = FALSE)



###########################################

# Cheching Benthic Master data

benthic_honduras<- read_csv(here("data", "raw", "BenthicMaster.csv")) %>% 
  filter(country == "Honduras")

# Duplicates: I think BenthicMaster is part of FF2.0! Check with George!

###########################################
# Check for duplicates
# 1. Get unique events from each source
events_by_source <- benthic_summary_combined %>%
  select(site, transect, sample_date, source) %>%
  distinct()

# 2. Spread into two sets
events_coralnet <- events_by_source %>% filter(source == "coralnet")
events_ff2 <- events_by_source %>% filter(source == "ff2")

# 3. Identify overlapping events
benthic_duplicated_events <- inner_join(
  events_coralnet,
  events_ff2,
  by = c("site", "transect", "sample_date")
)

# 4. View result
# nrow(benthic_duplicated_events)         # Number of duplicate entries
# benthic_duplicated_events               # Inspect duplicates

# There are 35 overlapping sampling events where both the coralnet dataset and the FF 2.0 dataset report benthic cover data for the same location, transect, and date. This indicates data duplication!! --> I will remove the coralnet ones since I think they are less reliable.


# Compare % cover of duplicates
# 1. Extract coralnet data for duplicated events
coralnet_dupes <- benthic_summary_combined %>%
  semi_join(benthic_duplicated_events %>% select(site, transect, sample_date),
            by = c("site", "transect", "sample_date")) %>%
  filter(source == "coralnet") %>%
  select(site, transect, sample_date, category, percent_cover_coralnet = percent_cover)

# 2. Extract ff2 data for duplicated events
ff2_dupes <- benthic_summary_combined %>%
  semi_join(benthic_duplicated_events %>% select(site, transect, sample_date),
            by = c("site", "transect", "sample_date")) %>%
  filter(source == "ff2") %>%
  select(site, transect, sample_date, category, percent_cover_ff2 = percent_cover)

# 3. Join both sources by site/transect/date/category
percent_cover_comparison <- full_join(
  coralnet_dupes,
  ff2_dupes,
  by = c("site", "transect", "sample_date", "category")
)

# 4. View result
# percent_cover_comparison


#Remove duplicates
benthic_summary_combined <- benthic_summary_combined %>%
  # Remove only the ff2 entries that match duplicated events
  anti_join(
    benthic_duplicated_events %>%
      select(site, transect, sample_date) %>%
      mutate(source = "ff2"),
    by = c("site", "transect", "sample_date", "source")
  )

#Available cleaned  data

# Fish data summary
fish_events_summary <- sample_events_df %>%
  mutate(year = year(sample_date)) %>%
  group_by(managed_access_areas, year) %>%
  summarise(fish_sample_events = n(), .groups = "drop") %>%
  arrange(managed_access_areas, year)

# Benthic data summary
benthic_events_summary <- benthic_summary_combined %>%
  select(site, transect, sample_date, managed_access_areas) %>%
  distinct() %>%
  mutate(year = year(sample_date)) %>%
  group_by(managed_access_areas, year) %>%
  summarise(benthic_sample_events = n(), .groups = "drop") %>%
  arrange(managed_access_areas, year)

# Join fish and benthic summaries by managed_access_areas and year and replace NAs with 0
sample_events_summary <- full_join(
  fish_events_summary,
  benthic_events_summary,
  by = c("managed_access_areas", "year")
) %>%
  mutate(
    fish_sample_events = replace_na(fish_sample_events, 0),
    benthic_sample_events = replace_na(benthic_sample_events, 0)
  ) %>%
  arrange(managed_access_areas, year)

# Sort the data
sample_events_summary <- sample_events_summary %>%
  arrange(managed_access_areas, year)

# Get the row numbers where group changes
break_rows <- which(
  sample_events_summary$managed_access_areas[-1] != sample_events_summary$managed_access_areas[-nrow(sample_events_summary)]
)

# Build the table
table_out <- sample_events_summary %>%
  kbl(
    caption = "Available data for Rare Honduras: sampling events from fish and benthic surveys",
    col.names = c("Managed Access Area", "Year", "Fish Sample Events", "Benthic Sample Events"),
    align = "lccc"
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Add horizontal lines after each group
for (row in break_rows) {
  table_out <- table_out %>% row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

table_out

Available data for Rare Honduras: sampling events from fish and benthic surveys

Managed Access Area	Year	Fish Sample Events	Benthic Sample Events
Belfate	2024	2	0
Guanaja	2019	13	0
Guanaja	2021	0	21
Guanaja	2022	6	12
Guanaja	2023	0	7
Iriona and Limon	2021	4	30
Puerto Cortes	2019	3	0
Puerto Cortes	2022	14	20
Roatan	2021	16	30
Roatan	2023	7	46
Santa Fe	2019	1	0
Santa Fe	2021	0	12
Santa Fe	2022	3	19
Santa Fe	2023	0	22
Trujillo	2019	3	0
Trujillo	2023	4	4
Utila	2023	6	15

---

Graphs

# Function to normalize management values
normalize_management_extended <- function(df) {
  df %>%
    mutate(management = case_when(
      management %in% c("Managed Access", "Managed Access Area", "Ley de Pesca y Acuicultura de Honduras") ~ "Managed Access",
      management %in% c("Reserve", "Zona de Recuperacion Pesquera", "Bajo Costa Azul", "La Venada", "Parque Nacional Marino Islas de la Bahia") ~ "Reserve",
      TRUE ~ management
    )) %>%
    filter(management != "Control") # Remove Control!!! (open access area)
}

# Apply to all three datasets
sample_events_df <- normalize_management_extended(sample_events_df)
sample_units_df <- normalize_management_extended(sample_units_df)
sample_observations_df <- normalize_management_extended(sample_observations_df)


# write.csv(sample_events_df, "sample_events_df.csv", row.names = FALSE)
# write.csv(sample_units_df, "sample_units_df.csv", row.names = FALSE)
# write.csv(sample_observations_df, "sample_observations_df.csv", row.names = FALSE)



# Check for average biomass based on the 3 df
# Aggregation from sample_events_df
events_avg <- sample_events_df %>%
  mutate(year = year(sample_date)) %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(biomass_kgha_avg_events = mean(biomass_kgha_avg, na.rm = TRUE), .groups = 'drop')

# Aggregation from sample_units_df
units_avg <- sample_units_df %>%
  mutate(year = year(sample_date)) %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(biomass_kgha_avg_units = mean(biomass_kgha, na.rm = TRUE), .groups = 'drop')

# Aggregation from sample_observations_df
# Step 1: Sum biomass per sample_unit_id
observations_unit <- sample_observations_df %>%
  group_by(sample_unit_id, sample_event_id, management, managed_access_areas, year = year(sample_date)) %>%
  summarise(biomass_unit = sum(biomass_kgha, na.rm = TRUE), .groups = 'drop')

# Step 2: Average these sums per sample_event_id
observations_event <- observations_unit %>%
  group_by(sample_event_id, management, managed_access_areas, year) %>%
  summarise(biomass_event = mean(biomass_unit, na.rm = TRUE), .groups = 'drop')

# Step 3: Annual average biomass by management and managed_access_areas
observations_avg <- observations_event %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(biomass_kgha_avg_observations = mean(biomass_event, na.rm = TRUE), .groups = 'drop')

# Combine all three aggregated datasets
biomass_comparison <- events_avg %>%
  left_join(units_avg, by = c("managed_access_areas", "management", "year")) %>%
  left_join(observations_avg, by = c("managed_access_areas", "management", "year"))

Fish Data

Fish Biomass Boxplots

1.1. Fish Biomass Boxplots Across Management Types and Years in Rare Honduras Sites

Interpretation: Each boxplot shows the distribution of fish biomass (t/km²) across different management types (Managed Access and Reserve) within a given Managed Access Area and year. The central horizontal line in each box represents the median biomass, while the top and bottom of the box indicate the interquartile range (IQR); i.e., the middle 50% of the data. The whiskers extend to the most extreme values within 1.5× the IQR, and any points beyond that are considered outliers. The individual black dots represent biomass values from individual transects, helping to visualize variability within each management type. These background points provide context on the number of observations and how spread out the biomass measurements are.

# Add fish plots:
# Aggregate observations_df correctly
observations_avg_df <- sample_observations_df %>%
  mutate(year = lubridate::year(sample_date)) %>%
  group_by(sample_unit_id, sample_event_id, managed_access_areas, management, year) %>%
  summarise(unit_biomass = sum(biomass_kgha, na.rm = TRUE), .groups = 'drop') %>%
  group_by(sample_event_id, managed_access_areas, management, year) %>%
  summarise(event_biomass = mean(unit_biomass, na.rm = TRUE), .groups = 'drop') %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(avg_biomass = mean(event_biomass, na.rm = TRUE), .groups = 'drop')


  
# Add year column
sample_events_df <- sample_events_df %>%
  mutate(year = lubridate::year(sample_date))

# Calculate mean and standard error per group
biomass_summary <- sample_events_df %>%
  group_by(managed_access_areas, year, management) %>%
  summarise(
    mean_biomass = mean(biomass_kgha_avg, na.rm = TRUE),
    se_biomass = sd(biomass_kgha_avg, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  )


# Prepare data: convert biomass from kg/ha to t/km²
units_df_prepped <- sample_units_df %>%
  mutate(
    year = lubridate::year(sample_date),
    biomass_t_km2 = biomass_kgha * 0.1  # conversion
  ) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019))

# Unique managed access areas
maas_units <- unique(units_df_prepped$managed_access_areas)

# Create plots
unit_plots <- list()

for (maa in maas_units) {
  df_subset <- units_df_prepped %>%
    filter(managed_access_areas == maa)
  
  # Determine upper limit for y-axis (98th percentile)
  upper_limit <- quantile(df_subset$biomass_t_km2, 0.98, na.rm = TRUE)

  p <- df_subset %>%
    ggplot(aes(x = management, y = biomass_t_km2, fill = management)) +
    geom_boxplot(alpha = 0.7, outlier.shape = NA) +
    geom_jitter(width = 0.2, alpha = 0.4, color = "black", size = 1) +
    facet_wrap(~ year) +
    coord_cartesian(ylim = c(0, upper_limit)) +
    scale_fill_manual(values = c("Reserve" = "#1B9E77", "Managed Access" = "#F8766D")) +
    labs(
      title = paste("Biomass (t/km²) by Management in", maa),
      x = "Management Type",
      y = "Biomass (t/km²)"
    ) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "none")
  
  unit_plots[[maa]] <- p
}

# Display plots
for (maa in maas_units) {
  print(unit_plots[[maa]])
}

# Median and IQT table
# Step 1: Compute summary table
biomass_summary_table <- units_df_prepped %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(
    median_biomass = median(biomass_kgha, na.rm = TRUE),
    lower_quartile = quantile(biomass_kgha, 0.25, na.rm = TRUE),
    upper_quartile = quantile(biomass_kgha, 0.75, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(managed_access_areas, management, year)

# Step 2: Round values and add t/km² columns
table_data <- biomass_summary_table %>%
  mutate(
    across(median_biomass:upper_quartile, ~ round(., 1)),
    median_t_km2 = round(median_biomass * 0.1, 2),
    lower_quartile_t_km2 = round(lower_quartile * 0.1, 2),
    upper_quartile_t_km2 = round(upper_quartile * 0.1, 2)
  )

# Step 3: Create table
table_out <- table_data %>%
  kbl(
    caption = "Biomass Summary by Managed Access Area, Management Type, and Year",
    col.names = c(
      "Managed Access Area", "Management Type", "Year",
      "Median Biomass (kg/ha)", "Lower Quartile (kg/ha)", "Upper Quartile (kg/ha)",
      "Median Biomass (t/km²)", "Lower Quartile (t/km²)", "Upper Quartile (t/km²)"
    )
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 4: Add horizontal lines after each province group
break_rows <- table_data %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 5: Display table
table_out

Biomass Summary by Managed Access Area, Management Type, and Year

Managed Access Area	Management Type	Year	Median Biomass (kg/ha)	Lower Quartile (kg/ha)	Upper Quartile (kg/ha)	Median Biomass (t/km²)	Lower Quartile (t/km²)	Upper Quartile (t/km²)
Belfate	Managed Access	2024	24.1	15.4	49.4	2.41	1.54	4.94
Guanaja	Managed Access	2019	90.5	36.1	170.2	9.05	3.61	17.02
Guanaja	Managed Access	2022	81.6	55.1	406.0	8.16	5.51	40.60
Guanaja	Reserve	2019	44.2	21.3	81.4	4.42	2.13	8.14
Guanaja	Reserve	2022	73.4	54.5	281.0	7.34	5.45	28.10
Iriona and Limon	Managed Access	2021	104.2	82.3	139.9	10.42	8.23	13.99
Puerto Cortes	Managed Access	2019	90.7	63.6	232.4	9.07	6.36	23.24
Puerto Cortes	Managed Access	2022	18.9	10.1	38.2	1.89	1.01	3.82
Puerto Cortes	Reserve	2022	25.9	8.4	64.0	2.59	0.84	6.40
Roatan	Managed Access	2021	157.6	61.4	265.7	15.76	6.14	26.57
Roatan	Reserve	2021	129.3	75.0	223.1	12.93	7.50	22.31
Roatan	Reserve	2023	162.2	71.3	226.3	16.22	7.13	22.63
Santa Fe	Managed Access	2022	290.7	155.4	425.9	29.07	15.54	42.59
Santa Fe	Reserve	2022	365.0	47.1	374.7	36.50	4.71	37.47
Trujillo	Managed Access	2019	63.6	55.5	118.8	6.36	5.55	11.88
Trujillo	Reserve	2023	47.0	35.3	176.2	4.70	3.53	17.62
Utila	Reserve	2023	91.0	53.0	164.0	9.10	5.30	16.40

1.2. Fish Biomass Boxplots per Managed Access Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Prepare data with unit conversion from kg/ha to t/km²
units_df_prepped <- sample_units_df %>%
  mutate(
    year = lubridate::year(sample_date),
    biomass_t_km2 = biomass_kgha * 0.1  # convert to t/km²
  ) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019))

# Step 2: Filter to only the most recent year per managed_access_areas
latest_years <- units_df_prepped %>%
  group_by(managed_access_areas) %>%
  summarise(latest_year = max(year, na.rm = TRUE), .groups = "drop")

units_df_latest <- units_df_prepped %>%
  inner_join(latest_years, by = "managed_access_areas") %>%
  filter(year == latest_year)

# Step 3: Get list of managed access areas (only those with latest year data)
maas_units <- unique(units_df_latest$managed_access_areas)

# Step 4: Create plots (latest year only)
unit_plots <- list()

for (maa in maas_units) {
  df_subset <- units_df_latest %>%
    filter(managed_access_areas == maa)

  upper_limit <- quantile(df_subset$biomass_t_km2, 0.98, na.rm = TRUE)

  p <- df_subset %>%
    ggplot(aes(x = factor(year), y = biomass_t_km2)) +
    geom_boxplot(alpha = 0.7, outlier.shape = NA, fill = "#6baed6") +
    geom_jitter(width = 0.2, alpha = 0.4, color = "black", size = 1) +
    coord_cartesian(ylim = c(0, upper_limit)) +
    labs(
      title = paste("Biomass (t/km²) in", maa),
      x = "Year",
      y = "Biomass (t/km²)"
    ) +
    theme_minimal(base_size = 12)

  unit_plots[[maa]] <- p
}

# Step 5: Display plots
for (maa in maas_units) {
  print(unit_plots[[maa]])
}

# Step 1: Compute the summary table without management
biomass_summary_table <- units_df_prepped %>%
  group_by(managed_access_areas, year) %>%
  summarise(
    median_biomass = median(biomass_kgha, na.rm = TRUE),
    lower_quartile = quantile(biomass_kgha, 0.25, na.rm = TRUE),
    upper_quartile = quantile(biomass_kgha, 0.75, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(managed_access_areas, year)

# Step 2: Round and add t/km² columns
table_data <- biomass_summary_table %>%
  mutate(
    across(median_biomass:upper_quartile, ~ round(., 1)),
    median_t_km2 = round(median_biomass * 0.1, 2),
    lower_quartile_t_km2 = round(lower_quartile * 0.1, 2),
    upper_quartile_t_km2 = round(upper_quartile * 0.1, 2)
  )

# Step 3: Create the kable table
table_out <- table_data %>%
  kbl(
    caption = "Biomass Summary by Managed Access Area and Year",
    col.names = c(
      "Managed Access Area", "Year",
      "Median Biomass (kg/ha)", "Lower Quartile (kg/ha)", "Upper Quartile (kg/ha)",
      "Median Biomass (t/km²)", "Lower Quartile (t/km²)", "Upper Quartile (t/km²)"
    )
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 4: Add horizontal lines after each Province
break_rows <- table_data %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 5: Display table
table_out

Biomass Summary by Managed Access Area and Year

Managed Access Area	Year	Median Biomass (kg/ha)	Lower Quartile (kg/ha)	Upper Quartile (kg/ha)	Median Biomass (t/km²)	Lower Quartile (t/km²)	Upper Quartile (t/km²)
Belfate	2024	24.1	15.4	49.4	2.41	1.54	4.94
Guanaja	2019	62.7	25.5	119.9	6.27	2.55	11.99
Guanaja	2022	74.6	55.1	364.4	7.46	5.51	36.44
Iriona and Limon	2021	104.2	82.3	139.9	10.42	8.23	13.99
Puerto Cortes	2019	90.7	63.6	232.4	9.07	6.36	23.24
Puerto Cortes	2022	21.9	8.4	60.4	2.19	0.84	6.04
Roatan	2021	147.1	71.4	256.7	14.71	7.14	25.67
Roatan	2023	162.2	71.3	226.3	16.22	7.13	22.63
Santa Fe	2022	365.0	42.3	408.9	36.50	4.23	40.89
Trujillo	2019	63.6	55.5	118.8	6.36	5.55	11.88
Trujillo	2023	47.0	35.3	176.2	4.70	3.53	17.62
Utila	2023	91.0	53.0	164.0	9.10	5.30	16.40

Mean Fish Biomass

2.1. Mean Fish Biomass with Standard Error by Management Type and Year in Rare Honduras Sites

Interpretation: Each bar in the chart represents the mean fish biomass (t/km²) for a given management type within a specific year and Managed Access Area. The height of the bar indicates the average biomass observed across sample events, while the black vertical line (error bar) represents the standard error of the mean, providing a sense of the variability and reliability of the estimate. Larger error bars suggest more variability or fewer observations. This visualization is helpful for comparing average biomass levels between management types over time and assessing the precision of these estimates.

# Step 1: Calculate mean biomass and standard error, convert to t/km²
units_summary_df <- sample_units_df %>%
  mutate(year = year(sample_date)) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019)) %>%
  group_by(managed_access_areas, year, management) %>%
  summarise(
    avg_biomass = mean(biomass_kgha, na.rm = TRUE) * 0.1,       # kg/ha → t/km²
    se_biomass = (sd(biomass_kgha, na.rm = TRUE) / sqrt(n())) * 0.1,
    .groups = "drop"
  )

# Step 2: Get unique managed access areas
maas_units <- unique(units_summary_df$managed_access_areas)

# Step 3: Create list of plots
unit_plots <- list()

for (maa in maas_units) {
  df_subset <- units_summary_df %>%
    filter(managed_access_areas == maa)

  upper_limit <- quantile(df_subset$avg_biomass + df_subset$se_biomass, 0.98, na.rm = TRUE)

  p <- ggplot(df_subset, aes(x = management, y = avg_biomass, fill = management)) +
    geom_col(alpha = 0.8, width = 0.7, position = position_dodge()) +
    geom_errorbar(
      aes(ymin = avg_biomass - se_biomass, ymax = avg_biomass + se_biomass),
      width = 0.2, position = position_dodge(0.7), color = "black"
    ) +
    facet_wrap(~ year) +
    coord_cartesian(ylim = c(0, upper_limit)) +
    scale_fill_manual(values = c("Reserve" = "#1B9E77", "Managed Access" = "#F8766D")) +
    labs(
      title = paste("Mean Biomass (t/km²) by Management in", maa),
      x = "Management Type",
      y = "Mean Biomass (t/km²)"
    ) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "none")

  unit_plots[[maa]] <- p
}

# Step 4: Print all plots
for (maa in maas_units) {
  print(unit_plots[[maa]])
}

units_summary_df <- sample_units_df %>%
  mutate(year = year(sample_date)) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019)) %>%
  group_by(managed_access_areas, year, management) %>%
  summarise(
    avg_biomass = mean(biomass_kgha, na.rm = TRUE),
    se_biomass = sd(biomass_kgha, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  )

# Add t/km² columns and round values
units_summary_table <- units_summary_df %>%
  mutate(
    avg_biomass = round(avg_biomass, 1),
    se_biomass = round(se_biomass, 1),
    avg_biomass_t_km2 = round(avg_biomass * 0.1, 2),
    se_biomass_t_km2 = round(se_biomass * 0.1, 2)
  )

# Select the correct columns
table_data <- units_summary_table %>%
  select(managed_access_areas, management, year,
         avg_biomass, se_biomass,
         avg_biomass_t_km2, se_biomass_t_km2)

# Format the table
table_out <- table_data %>%
  kbl(
    caption = "Mean Biomass and Standard Error by Managed Access Area, Management Type, and Year",
    col.names = c(
      "Managed Access Area", "Management Type", "Year", 
      "Mean Biomass (kg/ha)", "Standard Error (kg/ha)",
      "Mean Biomass (t/km²)", "Standard Error (t/km²)"
    ),
    align = "llcrrrr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Add horizontal lines after each province
break_rows <- table_data %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Show the table
table_out

Mean Biomass and Standard Error by Managed Access Area, Management Type, and Year

Managed Access Area	Management Type	Year	Mean Biomass (kg/ha)	Standard Error (kg/ha)	Mean Biomass (t/km²)	Standard Error (t/km²)
Belfate	Managed Access	2024	58.9	25.3	5.89	2.53
Guanaja	Managed Access	2019	113.0	13.9	11.30	1.39
Guanaja	Reserve	2019	65.9	9.3	6.59	0.93
Guanaja	Managed Access	2022	201.5	86.1	20.15	8.61
Guanaja	Reserve	2022	157.6	45.3	15.76	4.53
Iriona and Limon	Managed Access	2021	113.5	10.8	11.35	1.08
Puerto Cortes	Managed Access	2019	167.1	104.7	16.71	10.47
Puerto Cortes	Managed Access	2022	25.9	16.6	2.59	1.66
Puerto Cortes	Reserve	2022	50.0	11.0	5.00	1.10
Roatan	Managed Access	2021	250.6	80.2	25.06	8.02
Roatan	Reserve	2021	185.1	29.2	18.51	2.92
Roatan	Reserve	2023	167.4	18.3	16.74	1.83
Santa Fe	Managed Access	2022	290.7	270.4	29.07	27.04
Santa Fe	Reserve	2022	253.5	87.3	25.35	8.73
Trujillo	Managed Access	2019	95.0	39.8	9.50	3.98
Trujillo	Reserve	2023	143.1	45.2	14.31	4.52
Utila	Reserve	2023	142.8	32.6	14.28	3.26

2.2. Mean Fish Biomass with Standard Error per Managed Access Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Create summary table with all years
units_summary_df <- sample_units_df %>%
  mutate(year = year(sample_date)) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019)) %>%
  group_by(managed_access_areas, year) %>%
  summarise(
    avg_biomass = mean(biomass_kgha, na.rm = TRUE) * 0.1,  # t/km²
    se_biomass = (sd(biomass_kgha, na.rm = TRUE) / sqrt(n())) * 0.1,
    .groups = "drop"
  )

# Step 2: Use this to create the table (all years)
units_summary_table <- units_summary_df %>%
  mutate(
    avg_biomass = round(avg_biomass, 1),
    se_biomass = round(se_biomass, 1),
    avg_biomass_t_km2 = round(avg_biomass, 2),  # Already in t/km²
    se_biomass_t_km2 = round(se_biomass, 2)
  )

# Step 3: Create version with only latest year (for plotting)
latest_years <- units_summary_df %>%
  group_by(managed_access_areas) %>%
  summarise(latest_year = max(year, na.rm = TRUE), .groups = "drop")

units_summary_df_latest <- units_summary_df %>%
  inner_join(latest_years, by = "managed_access_areas") %>%
  filter(year == latest_year)

# Step 4: Get list of MAAs
maas_units <- unique(units_summary_df_latest$managed_access_areas)

# Step 5: Generate plots
unit_plots <- list()

for (maa in maas_units) {
  df_subset <- units_summary_df_latest %>%
    filter(managed_access_areas == maa)

  upper_limit <- quantile(df_subset$avg_biomass + df_subset$se_biomass, 0.98, na.rm = TRUE)

  p <- ggplot(df_subset, aes(x = factor(year), y = avg_biomass)) +
    geom_col(fill = "#6baed6", alpha = 0.8, width = 0.7) +
    geom_errorbar(
      aes(ymin = avg_biomass - se_biomass, ymax = avg_biomass + se_biomass),
      width = 0.2, color = "black"
    ) +
    coord_cartesian(ylim = c(0, upper_limit)) +
    labs(
      title = paste("Mean Biomass (t/km²) in", maa),
      x = "Year",
      y = "Mean Biomass (t/km²)"
    ) +
    theme_minimal(base_size = 12)

  unit_plots[[maa]] <- p
}

# Step 6: Show plots
for (maa in maas_units) {
  print(unit_plots[[maa]])
}

# Step 1: Add properly rounded values for display in both units
units_summary_table <- units_summary_df %>%
  mutate(
    avg_biomass_t_km2 = round(avg_biomass, 2),
    se_biomass_t_km2 = round(se_biomass, 2),
    avg_biomass_kg_ha = round(avg_biomass * 100, 1),  # display in kg/ha
    se_biomass_kg_ha = round(se_biomass * 100, 1)
  )

# Step 2: Select correct columns (no management)
table_data <- units_summary_table %>%
  select(managed_access_areas, year,
         avg_biomass_kg_ha, se_biomass_kg_ha,
         avg_biomass_t_km2, se_biomass_t_km2)

# Step 3: Format the table
table_out <- table_data %>%
  kbl(
    caption = "Mean Biomass and Standard Error by Managed Access Area and Year",
    col.names = c(
      "Managed Access Area", "Year",
      "Mean Biomass (kg/ha)", "Standard Error (kg/ha)",
      "Mean Biomass (t/km²)", "Standard Error (t/km²)"
    ),
    align = "llrrrr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 4: Add horizontal separators after each province
break_rows <- table_data %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 5: Show table
table_out

Mean Biomass and Standard Error by Managed Access Area and Year

Managed Access Area	Year	Mean Biomass (kg/ha)	Standard Error (kg/ha)	Mean Biomass (t/km²)	Standard Error (t/km²)
Belfate	2024	589.1	252.9	5.89	2.53
Guanaja	2019	886.6	85.0	8.87	0.85
Guanaja	2022	1704.7	395.7	17.05	3.96
Iriona and Limon	2021	1135.1	107.5	11.35	1.08
Puerto Cortes	2019	1671.4	1047.0	16.71	10.47
Puerto Cortes	2022	480.2	101.8	4.80	1.02
Roatan	2021	2119.9	370.4	21.20	3.70
Roatan	2023	1674.3	183.3	16.74	1.83
Santa Fe	2022	2641.3	846.0	26.41	8.46
Trujillo	2019	950.2	397.8	9.50	3.98
Trujillo	2023	1430.7	452.1	14.31	4.52
Utila	2023	1427.5	326.3	14.28	3.26

Mean Fish Biomass by Trophic Group

3.1. Mean Fish Biomass by Trophic Group, Management Type, and Year in Rare Honduras Sites

Interpretation: This bar chart displays the mean fish biomass (t/km²) for different trophic groups across management types and years. Each bar represents the average biomass within a specific trophic group and management type, while the black vertical lines indicate the standard error, reflecting variability across sampling events.

# Step 1: Simplify trophic groups and add year
sample_observations_df <- sample_observations_df %>%
  mutate(
    year = lubridate::year(sample_date),
    trophic_group = case_when(
      trophic_group %in% c("herbivore-macroalgae", "herbivore-detritivore") ~ "herbivore",
      trophic_group %in% c("invertivore-mobile", "invertivore-sessile") ~ "invertivore",
      TRUE ~ trophic_group
    )
  )

# Step 2: Sum biomass per unit (still in kg/ha)
obs_unit_biomass <- sample_observations_df %>%
  group_by(sample_unit_id, sample_event_id, managed_access_areas, management, trophic_group, year) %>%
  summarise(unit_biomass = sum(biomass_kgha, na.rm = TRUE), .groups = "drop")

# Step 3: Average per event
obs_event_biomass <- obs_unit_biomass %>%
  group_by(sample_event_id, managed_access_areas, management, trophic_group, year) %>%
  summarise(event_biomass = mean(unit_biomass, na.rm = TRUE), .groups = "drop")

# Step 4: Calculate mean and SE in t/km²
corrected_trophic_summary <- obs_event_biomass %>%
  group_by(managed_access_areas, year, management, trophic_group) %>%
  summarise(
    mean_biomass = mean(event_biomass, na.rm = TRUE) * 0.1,  # convert to t/km²
    se_biomass = (sd(event_biomass, na.rm = TRUE) / sqrt(n())) * 0.1,
    .groups = "drop"
  ) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019))

# Step 5: Plot generation
maas <- unique(corrected_trophic_summary$managed_access_areas)
trophic_plots <- list()

trophic_levels <- c("herbivore", "invertivore", "omnivore", "piscivore", "planktivore")

for (maa in maas) {
  df_subset <- corrected_trophic_summary %>%
    filter(managed_access_areas == maa) %>%
    mutate(trophic_group = factor(trophic_group, levels = trophic_levels))

  p <- ggplot(df_subset, aes(x = trophic_group, y = mean_biomass, fill = management)) +
    geom_col(position = position_dodge(width = 0.9), width = 0.8) +
    geom_errorbar(
      aes(ymin = mean_biomass - se_biomass, ymax = mean_biomass + se_biomass),
      position = position_dodge(width = 0.9), width = 0.2
    ) +
    facet_wrap(~ year) +
    scale_fill_manual(values = c("Reserve" = "#1B9E77", "Managed Access" = "#F8766D")) +
    labs(
      title = paste("Mean Biomass (t/km²) by Trophic Group in", maa),
      x = "Trophic Group",
      y = "Mean Biomass (t/km²)",
      fill = "Management"
    ) +
    theme_minimal(base_size = 13) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  trophic_plots[[maa]] <- p
}

# Step 6: Display plots
for (maa in maas) {
  print(trophic_plots[[maa]])
}

3.2. Mean Fish Biomass by Trophic Group, per Managed Acces Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Prepare trophic groups and extract year
sample_observations_df <- sample_observations_df %>%
  mutate(
    year = lubridate::year(sample_date),
    trophic_group = case_when(
      trophic_group %in% c("herbivore-macroalgae", "herbivore-detritivore") ~ "herbivore",
      trophic_group %in% c("invertivore-mobile", "invertivore-sessile") ~ "invertivore",
      TRUE ~ trophic_group
    )
  )

# Step 2: Sum biomass per unit
obs_unit_biomass <- sample_observations_df %>%
  group_by(sample_unit_id, sample_event_id, managed_access_areas, trophic_group, year) %>%
  summarise(unit_biomass = sum(biomass_kgha, na.rm = TRUE), .groups = "drop")

# Step 3: Average per event
obs_event_biomass <- obs_unit_biomass %>%
  group_by(sample_event_id, managed_access_areas, trophic_group, year) %>%
  summarise(event_biomass = mean(unit_biomass, na.rm = TRUE), .groups = "drop")

# Step 4: Calculate mean and SE in t/km²
corrected_trophic_summary_all_years <- obs_event_biomass %>%
  group_by(managed_access_areas, year, trophic_group) %>%
  summarise(
    mean_biomass = mean(event_biomass, na.rm = TRUE) * 0.1,  # kg/ha → t/km²
    se_biomass = (sd(event_biomass, na.rm = TRUE) / sqrt(n())) * 0.1,
    .groups = "drop"
  ) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019))

# Step 5: Filter to only the latest year per MAA
latest_years <- corrected_trophic_summary_all_years %>%
  group_by(managed_access_areas) %>%
  summarise(latest_year = max(year), .groups = "drop")

corrected_trophic_summary <- corrected_trophic_summary_all_years %>%
  inner_join(latest_years, by = c("managed_access_areas", "year" = "latest_year"))

# Step 6: Plotting (no faceting now, just one year per MAA)
maas <- unique(corrected_trophic_summary$managed_access_areas)
trophic_plots <- list()

trophic_levels <- c("herbivore", "invertivore", "omnivore", "piscivore", "planktivore")

for (maa in maas) {
  df_subset <- corrected_trophic_summary %>%
    filter(managed_access_areas == maa) %>%
    mutate(trophic_group = factor(trophic_group, levels = trophic_levels))

  upper_limit <- quantile(df_subset$mean_biomass + df_subset$se_biomass, 0.98, na.rm = TRUE)

  p <- ggplot(df_subset, aes(x = trophic_group, y = mean_biomass)) +
    geom_col(fill = "#6baed6", alpha = 0.8, width = 0.8) +
    geom_errorbar(
      aes(ymin = mean_biomass - se_biomass, ymax = mean_biomass + se_biomass),
      width = 0.2
    ) +
    labs(
      title = paste("Mean Biomass (t/km²) by Trophic Group in", maa, "-", unique(df_subset$year)),
      x = "Trophic Group",
      y = "Mean Biomass (t/km²)"
    ) +
    coord_cartesian(ylim = c(0, upper_limit)) +
    theme_minimal(base_size = 13) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  trophic_plots[[maa]] <- p
}

# Step 7: Print all plots
for (maa in maas) {
  print(trophic_plots[[maa]])
}

Mean Herbivore Fish Biomass

4.1. Mean Total and Herbivore Fish Biomass by Management Type and Year in Rare Honduras Sites

Interpretation: This chart presents the mean total fish biomass (t/km²) per management type and year, along with a highlighted overlay of herbivore biomass. The solid bars represent the average biomass of herbivorous fish, while the lighter shaded bars extending above each bar indicate the average total biomass. The black error bars represent the standard error of the average herbivorous fish biomass, showing variability across observations. This visualization allows for a direct comparison between management types and years, and helps assess the relative contribution of herbivores to the total biomass, a useful indicator of reef health and fishing pressure.

# Prepare simplified trophic groups and standardize fields
sample_observations_df <- sample_observations_df %>%
  mutate(
    trophic_group_simplified = case_when(
      trophic_group %in% c("herbivore", "herbivore-macroalgae", "herbivore-detritivore") ~ "Herbivore",
      TRUE ~ "Other"
    ),
    year = lubridate::year(sample_date),
    management = case_when(
      management %in% c("Managed Access", "Managed Access Area", "Ley de Pesca y Acuicultura de Honduras") ~ "Managed Access",
      management %in% c("Reserve", "Zona de Recuperacion Pesquera") ~ "Reserve",
      TRUE ~ management
    ),
    trophic_group_simplified = factor(trophic_group_simplified, levels = c("Herbivore", "Other"))
  )

# Step 1: Sum biomass per unit
herb_unit_biomass <- sample_observations_df %>%
  group_by(sample_unit_id, sample_event_id, managed_access_areas, management, trophic_group_simplified, year) %>%
  summarise(unit_biomass = sum(biomass_kgha, na.rm = TRUE), .groups = "drop")

# Step 2: Average per event
herb_event_biomass <- herb_unit_biomass %>%
  group_by(sample_event_id, managed_access_areas, management, trophic_group_simplified, year) %>%
  summarise(event_biomass = mean(unit_biomass, na.rm = TRUE), .groups = "drop")

# Step 3: Compute mean and SE in t/km²
herb_summary_avg <- herb_event_biomass %>%
  group_by(managed_access_areas, management, year, trophic_group_simplified) %>%
  summarise(
    avg_biomass = mean(event_biomass, na.rm = TRUE) * 0.1,
    se_biomass = (sd(event_biomass, na.rm = TRUE) / sqrt(n())) * 0.1,
    .groups = "drop"
  ) %>%
  mutate(se_biomass = ifelse(is.na(se_biomass), 0, se_biomass)) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019))

# Plot: Herbivore (with SE) vs. Other (transparent)
herb_stacked_plots <- list()

for (maa in unique(herb_summary_avg$managed_access_areas)) {
  df_subset <- herb_summary_avg %>%
    filter(managed_access_areas == maa) %>%
    mutate(
      alpha_val = ifelse(trophic_group_simplified == "Other", 0.4, 1)
    )

  df_herbivore <- df_subset %>% filter(trophic_group_simplified == "Herbivore")
  dummy_legend <- data.frame(x = 1, y = 1, label = "Herbivore Biomass")

  p <- ggplot(df_subset, aes(x = management, y = avg_biomass, fill = management, alpha = alpha_val)) +
    geom_col(position = "stack", width = 0.8, color = NA) +

    geom_errorbar(
  data = df_herbivore,
  aes(x = management, ymin = avg_biomass - se_biomass, ymax = avg_biomass + se_biomass),
  width = 0.2,
  position = position_dodge(width = 0.8),
  color = "black",
  inherit.aes = FALSE
) +

    geom_col(
  data = df_herbivore,
  aes(x = management, y = avg_biomass),
  fill = NA,
  color = "black",
  position = "stack",
  width = 0.8,
  linewidth = 0.5,
  inherit.aes = FALSE
) +

    ggnewscale::new_scale("shape") +
    geom_point(
      data = dummy_legend,
      mapping = aes(x = x, y = y, shape = label),
      size = NA,
      color = "black",
      fill = "white",
      stroke = 0.7,
      inherit.aes = FALSE,
      show.legend = TRUE
    ) +
    scale_shape_manual(
      name = "",
      values = c("Herbivore Biomass" = 22),
      guide = guide_legend(override.aes = list(
        shape = 22,
        fill = "white",
        color = "black",
        size = 4
      ))
    ) +

    facet_wrap(~year) +
    scale_fill_manual(
      values = c("Reserve" = "#1B9E77", "Managed Access" = "#F8766D"),
      guide = "none"
    ) +
    scale_alpha_identity() +
    labs(
      title = paste("Mean Biomass by Management in", maa),
      x = "Management Type",
      y = "Mean Biomass (t/km²)"
    ) +
    theme_minimal(base_size = 12) +
    theme(axis.text.x = element_text(angle = 0))

  herb_stacked_plots[[maa]] <- p
}

# Show plots
for (maa in unique(herb_summary_avg$managed_access_areas)) {
  print(herb_stacked_plots[[maa]])
}

# Step 1: Filter herbivore and compute both units correctly
herbivore_table <- herb_summary_avg %>%
  filter(trophic_group_simplified == "Herbivore") %>%
  mutate(
    avg_biomass_t_km2 = round(avg_biomass, 2),                   # already in t/km²
    se_biomass_t_km2 = round(se_biomass, 2),
    avg_biomass_kg_ha = round(avg_biomass * 100, 1),             # convert to kg/ha
    se_biomass_kg_ha = round(se_biomass * 100, 1)
  )

# Step 2: Format table with both units
table_out <- herbivore_table %>%
  select(
    managed_access_areas, management, year,
    avg_biomass_kg_ha, se_biomass_kg_ha,
    avg_biomass_t_km2, se_biomass_t_km2
  ) %>%
  kbl(
    caption = "Mean Herbivore Biomass and Standard Error by Managed Access Area, Management Type, and Year",
    col.names = c(
      "Managed Access Area", "Management Type", "Year",
      "Mean Biomass (kg/ha)", "Standard Error (kg/ha)",
      "Mean Biomass (t/km²)", "Standard Error (t/km²)"
    ),
    align = "lllrrrr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 3: Add horizontal lines by MAA
break_rows <- herbivore_table %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 4: Show table
table_out

Mean Herbivore Biomass and Standard Error by Managed Access Area, Management Type, and Year

Managed Access Area	Management Type	Year	Mean Biomass (kg/ha)	Standard Error (kg/ha)	Mean Biomass (t/km²)	Standard Error (t/km²)
Belfate	Managed Access	2024	26.2	0.0	0.26	0.00
Guanaja	Managed Access	2019	838.3	181.4	8.38	1.81
Guanaja	Managed Access	2022	925.1	375.6	9.25	3.76
Guanaja	Reserve	2019	690.6	310.3	6.91	3.10
Guanaja	Reserve	2022	462.3	143.4	4.62	1.43
Iriona and Limon	Managed Access	2021	181.8	43.6	1.82	0.44
Puerto Cortes	Managed Access	2019	566.4	205.0	5.66	2.05
Puerto Cortes	Managed Access	2022	216.6	0.0	2.17	0.00
Puerto Cortes	Reserve	2022	302.2	93.3	3.02	0.93
Roatan	Managed Access	2021	1406.6	615.4	14.07	6.15
Roatan	Reserve	2021	1233.3	272.2	12.33	2.72
Roatan	Reserve	2023	1111.3	144.2	11.11	1.44
Santa Fe	Managed Access	2022	318.4	265.8	3.18	2.66
Santa Fe	Reserve	2022	1235.9	0.0	12.36	0.00
Trujillo	Managed Access	2019	510.8	431.8	5.11	4.32
Trujillo	Reserve	2023	713.1	240.0	7.13	2.40
Utila	Reserve	2023	519.0	21.8	5.19	0.22

4.2. Mean Total and Herbivore Fish Biomass per Managed Access Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Simplify trophic groups and extract year
sample_observations_df <- sample_observations_df %>%
  mutate(
    trophic_group_simplified = case_when(
      trophic_group %in% c("herbivore", "herbivore-macroalgae", "herbivore-detritivore") ~ "Herbivore",
      TRUE ~ "Other"
    ),
    year = lubridate::year(sample_date),
    trophic_group_simplified = factor(trophic_group_simplified, levels = c("Herbivore", "Other"))
  )

# Step 2: Sum biomass per unit
herb_unit_biomass <- sample_observations_df %>%
  group_by(sample_unit_id, sample_event_id, managed_access_areas, trophic_group_simplified, year) %>%
  summarise(unit_biomass = sum(biomass_kgha, na.rm = TRUE), .groups = "drop")

# Step 3: Average per event
herb_event_biomass <- herb_unit_biomass %>%
  group_by(sample_event_id, managed_access_areas, trophic_group_simplified, year) %>%
  summarise(event_biomass = mean(unit_biomass, na.rm = TRUE), .groups = "drop")

# Step 4: Mean & SE per year/MAA/trophic group — now in t/km²
herb_summary_avg <- herb_event_biomass %>%
  group_by(managed_access_areas, year, trophic_group_simplified) %>%
  summarise(
    avg_biomass = mean(event_biomass, na.rm = TRUE) * 0.1,
    se_biomass = (sd(event_biomass, na.rm = TRUE) / sqrt(n())) * 0.1,
    .groups = "drop"
  ) %>%
  mutate(se_biomass = ifelse(is.na(se_biomass), 0, se_biomass)) %>%
  filter(!(managed_access_areas == "Santa Fe" & year == 2019))

# Step 5 (updated): Create a filtered version for plotting only the latest year
latest_years <- herb_summary_avg %>%
  group_by(managed_access_areas) %>%
  summarise(latest_year = max(year, na.rm = TRUE), .groups = "drop")

herb_summary_latest <- herb_summary_avg %>%
  inner_join(latest_years, by = c("managed_access_areas", "year" = "latest_year"))

# Step 6 (updated): Plot — Stacked bars for latest year only
herb_stacked_plots <- list()

for (maa in unique(herb_summary_latest$managed_access_areas)) {
  df_subset <- herb_summary_latest %>%
    filter(managed_access_areas == maa) %>%
    mutate(alpha_val = ifelse(trophic_group_simplified == "Other", 0.4, 1))

  df_herbivore <- df_subset %>%
    filter(trophic_group_simplified == "Herbivore")

  dummy_legend <- data.frame(x = 1, y = 1, label = "Herbivore Biomass")

  p <- ggplot(df_subset, aes(x = factor(year), y = avg_biomass, alpha = alpha_val)) +
    geom_col(position = "stack", width = 0.8, fill = "#6baed6", color = NA) +

    geom_errorbar(
      data = df_herbivore,
      aes(x = factor(year), ymin = avg_biomass - se_biomass, ymax = avg_biomass + se_biomass),
      width = 0.2,
      color = "black",
      inherit.aes = FALSE
    ) +

    geom_col(
      data = df_herbivore,
      aes(x = factor(year), y = avg_biomass),
      fill = NA,
      color = "black",
      position = "stack",
      width = 0.8,
      linewidth = 0.5,
      inherit.aes = FALSE
    ) +

    ggnewscale::new_scale("shape") +
    geom_point(
      data = dummy_legend,
      mapping = aes(x = x, y = y, shape = label),
      size = NA,
      color = "black",
      fill = "white",
      stroke = 0.7,
      inherit.aes = FALSE,
      show.legend = TRUE
    ) +
    scale_shape_manual(
      name = "",
      values = c("Herbivore Biomass" = 22),
      guide = guide_legend(override.aes = list(
        shape = 22,
        fill = "white",
        color = "black",
        size = 4
      ))
    ) +

    scale_alpha_identity() +
    labs(
      title = paste("Mean Biomass by Trophic Group in", maa),
      x = "Year",
      y = "Mean Biomass (t/km²)"
    ) +
    theme_minimal(base_size = 12) +
    theme(axis.text.x = element_text(angle = 0))

  herb_stacked_plots[[maa]] <- p
}

# Step 7: Display all plots
for (maa in unique(herb_summary_latest$managed_access_areas)) {
  print(herb_stacked_plots[[maa]])
}

# Step 1: Filter Herbivore rows and create both kg/ha and t/km² values properly
herbivore_table <- herb_summary_avg %>%
  filter(trophic_group_simplified == "Herbivore") %>%
  mutate(
    avg_biomass_t_km2 = round(avg_biomass, 2),
    se_biomass_t_km2 = round(se_biomass, 2),
    avg_biomass_kg_ha = round(avg_biomass * 100, 1),
    se_biomass_kg_ha = round(se_biomass * 100, 1)
  )

# Step 2: Create formatted table
table_data <- herbivore_table %>%
  select(managed_access_areas, year,
         avg_biomass_kg_ha, se_biomass_kg_ha,
         avg_biomass_t_km2, se_biomass_t_km2)

table_out <- table_data %>%
  kbl(
    caption = "Mean Herbivore Biomass and Standard Error by Managed Access Area and Year",
    col.names = c(
      "Managed Access Area", "Year",
      "Mean Herbivore Biomass (kg/ha)", "Standard Error (kg/ha)",
      "Mean Herbivore Biomass (t/km²)", "Standard Error (t/km²)"
    ),
    align = "llrrrr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 3: Add horizontal lines after each Province
break_rows <- table_data %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 4: Show the table
table_out

Mean Herbivore Biomass and Standard Error by Managed Access Area and Year

Managed Access Area	Year	Mean Herbivore Biomass (kg/ha)	Standard Error (kg/ha)	Mean Herbivore Biomass (t/km²)	Standard Error (t/km²)
Belfate	2024	26.2	0.0	0.26	0.00
Guanaja	2019	758.8	180.8	7.59	1.81
Guanaja	2022	616.6	164.8	6.17	1.65
Iriona and Limon	2021	181.8	43.6	1.82	0.44
Puerto Cortes	2019	566.4	205.0	5.66	2.05
Puerto Cortes	2022	295.6	86.1	2.96	0.86
Roatan	2021	1309.1	298.3	13.09	2.98
Roatan	2023	1111.3	144.2	11.11	1.44
Santa Fe	2022	624.2	342.2	6.24	3.42
Trujillo	2019	510.8	431.8	5.11	4.32
Trujillo	2023	713.1	240.0	7.13	2.40
Utila	2023	519.0	21.8	5.19	0.22

---

Benthic Data

Mean Benthic Cover

1.1. Mean Benthic Cover Composition by Management Type and Year in Rare Honduras Sites

Interpretation: This stacked bar chart shows the mean percent cover of various benthic groups across Managed Access and Reserve sites for multiple years and across all Managed Access Areas. Each bar represents the total benthic community composition for a given year and management type, with each colored segment indicating the average contribution of a specific benthic group.

# Normalize and summarize per transect before aggregating
benthic_df_prepped <- benthic_summary_combined %>%
  mutate(year = year(sample_date)) %>%
  group_by(year, management, managed_access_areas, category) %>%
  summarise(
    mean_raw_cover = mean(percent_cover, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(year, management, managed_access_areas) %>%
  mutate(
    total_cover = sum(mean_raw_cover),
    normalized_cover = 100 * mean_raw_cover / total_cover
  ) %>%
  ungroup()

# List of plots
maas <- unique(benthic_df_prepped$managed_access_areas)
benthic_composition_stacked <- list()

for (maa in maas) {
  df_subset <- benthic_df_prepped %>%
    filter(managed_access_areas == maa)
  
  benthic_colors <- c(
    "Hard coral" = "#C77CFF",
    "Soft coral" = "#F8766D",
    "Macroalgae" = "#66A61E",
    "Turf algae" = "#1B9E77",
    "Other invertebrates" = "#E6AB02",
    "Rubble" = "#A6761D",
    "Sand" = "#FFD700",
    "Bare substrate" = "#666666",
    "Other" = "#999999"
  )
  
  p <- ggplot(df_subset, aes(x = management, y = normalized_cover, fill = category)) +
    geom_bar(stat = "identity", position = "stack", width = 0.7, alpha = 0.8) +
    facet_wrap(~ year) +
    labs(
      title = paste("Mean Benthic Cover Composition in", maa),
      x = "Management Type",
      y = "Mean Percent Cover",
      fill = "Benthic Group"
    ) +
    scale_fill_manual(values = benthic_colors) +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "bottom",
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9),
      legend.key.size = unit(0.5, "lines")  # controls box size
    )

  benthic_composition_stacked[[maa]] <- p
}


for (maa in maas) {
  print(benthic_composition_stacked[[maa]])
}

# Step 1: Prepare and round the data
benthic_table <- benthic_df_prepped %>%
  group_by(managed_access_areas, management, year, category) %>%
  summarise(normalized_cover = mean(normalized_cover, na.rm = TRUE), .groups = "drop") %>%
  mutate(
    normalized_cover = round(normalized_cover, 0)
  ) %>%
  select(managed_access_areas, management, year, category, normalized_cover)

# Step 2: Format the table with kableExtra
table_out <- benthic_table %>%
  kbl(
    caption = "Benthic Cover (%) by Managed Access Area, Management Type, and Year",
    col.names = c("Managed Access Area", "Management Type", "Year", "Benthic Category", "Percent Cover"),
    align = "lllcr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 3: Add horizontal lines after each Province/Management/Year group
break_rows <- benthic_table %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 4: Show the final table
table_out

Benthic Cover (%) by Managed Access Area, Management Type, and Year

Managed Access Area	Management Type	Year	Benthic Category	Percent Cover
Guanaja	Managed Access	2021	Hard coral	6
Guanaja	Managed Access	2021	Macroalgae	17
Guanaja	Managed Access	2021	Other	54
Guanaja	Managed Access	2021	Other invertebrates	6
Guanaja	Managed Access	2021	Sand	17
Guanaja	Managed Access	2022	Hard coral	4
Guanaja	Managed Access	2022	Macroalgae	33
Guanaja	Managed Access	2022	Other	15
Guanaja	Managed Access	2022	Other invertebrates	2
Guanaja	Managed Access	2022	Rubble	5
Guanaja	Managed Access	2022	Sand	11
Guanaja	Managed Access	2022	Soft coral	11
Guanaja	Managed Access	2022	Turf algae	20
Guanaja	Reserve	2021	Hard coral	1
Guanaja	Reserve	2021	Macroalgae	6
Guanaja	Reserve	2021	Other	86
Guanaja	Reserve	2021	Other invertebrates	2
Guanaja	Reserve	2021	Sand	6
Guanaja	Reserve	2022	Hard coral	4
Guanaja	Reserve	2022	Macroalgae	33
Guanaja	Reserve	2022	Other	19
Guanaja	Reserve	2022	Other invertebrates	2
Guanaja	Reserve	2022	Rubble	2
Guanaja	Reserve	2022	Sand	6
Guanaja	Reserve	2022	Soft coral	18
Guanaja	Reserve	2022	Turf algae	16
Guanaja	Reserve	2023	Hard coral	9
Guanaja	Reserve	2023	Macroalgae	39
Guanaja	Reserve	2023	Other	11
Guanaja	Reserve	2023	Other invertebrates	1
Guanaja	Reserve	2023	Rubble	1
Guanaja	Reserve	2023	Sand	9
Guanaja	Reserve	2023	Soft coral	10
Guanaja	Reserve	2023	Turf algae	20
Iriona and Limon	Managed Access	2021	Bare substrate	1
Iriona and Limon	Managed Access	2021	Hard coral	8
Iriona and Limon	Managed Access	2021	Macroalgae	20
Iriona and Limon	Managed Access	2021	Other	3
Iriona and Limon	Managed Access	2021	Other invertebrates	9
Iriona and Limon	Managed Access	2021	Sand	48
Iriona and Limon	Managed Access	2021	Soft coral	4
Iriona and Limon	Managed Access	2021	Turf algae	8
Iriona and Limon	Reserve	2021	Hard coral	10
Iriona and Limon	Reserve	2021	Macroalgae	65
Iriona and Limon	Reserve	2021	Other	2
Iriona and Limon	Reserve	2021	Other invertebrates	6
Iriona and Limon	Reserve	2021	Rubble	1
Iriona and Limon	Reserve	2021	Sand	7
Iriona and Limon	Reserve	2021	Soft coral	4
Iriona and Limon	Reserve	2021	Turf algae	6
Puerto Cortes	Reserve	2022	Hard coral	3
Puerto Cortes	Reserve	2022	Macroalgae	24
Puerto Cortes	Reserve	2022	Other	4
Puerto Cortes	Reserve	2022	Other invertebrates	3
Puerto Cortes	Reserve	2022	Sand	8
Puerto Cortes	Reserve	2022	Soft coral	6
Puerto Cortes	Reserve	2022	Turf algae	52
Roatan	Managed Access	2021	Bare substrate	8
Roatan	Managed Access	2021	Hard coral	15
Roatan	Managed Access	2021	Macroalgae	12
Roatan	Managed Access	2021	Other	4
Roatan	Managed Access	2021	Other invertebrates	5
Roatan	Managed Access	2021	Rubble	2
Roatan	Managed Access	2021	Sand	8
Roatan	Managed Access	2021	Soft coral	10
Roatan	Managed Access	2021	Turf algae	35
Roatan	Managed Access	2023	Bare substrate	1
Roatan	Managed Access	2023	Hard coral	5
Roatan	Managed Access	2023	Macroalgae	20
Roatan	Managed Access	2023	Other	40
Roatan	Managed Access	2023	Other invertebrates	2
Roatan	Managed Access	2023	Rubble	3
Roatan	Managed Access	2023	Sand	5
Roatan	Managed Access	2023	Soft coral	7
Roatan	Managed Access	2023	Turf algae	16
Roatan	Reserve	2021	Bare substrate	4
Roatan	Reserve	2021	Hard coral	16
Roatan	Reserve	2021	Macroalgae	10
Roatan	Reserve	2021	Other	7
Roatan	Reserve	2021	Other invertebrates	3
Roatan	Reserve	2021	Rubble	2
Roatan	Reserve	2021	Sand	10
Roatan	Reserve	2021	Soft coral	15
Roatan	Reserve	2021	Turf algae	32
Santa Fe	Managed Access	2022	Hard coral	9
Santa Fe	Managed Access	2022	Macroalgae	3
Santa Fe	Managed Access	2022	Other invertebrates	81
Santa Fe	Managed Access	2022	Sand	3
Santa Fe	Managed Access	2022	Soft coral	2
Santa Fe	Managed Access	2022	Turf algae	1
Santa Fe	Managed Access	2023	Hard coral	8
Santa Fe	Managed Access	2023	Macroalgae	23
Santa Fe	Managed Access	2023	Other	5
Santa Fe	Managed Access	2023	Other invertebrates	42
Santa Fe	Managed Access	2023	Rubble	0
Santa Fe	Managed Access	2023	Sand	9
Santa Fe	Managed Access	2023	Soft coral	2
Santa Fe	Managed Access	2023	Turf algae	10
Santa Fe	Reserve	2021	Bare substrate	1
Santa Fe	Reserve	2021	Hard coral	8
Santa Fe	Reserve	2021	Macroalgae	40
Santa Fe	Reserve	2021	Other	3
Santa Fe	Reserve	2021	Other invertebrates	14
Santa Fe	Reserve	2021	Rubble	0
Santa Fe	Reserve	2021	Sand	9
Santa Fe	Reserve	2021	Soft coral	5
Santa Fe	Reserve	2021	Turf algae	20
Santa Fe	Reserve	2022	Hard coral	12
Santa Fe	Reserve	2022	Macroalgae	37
Santa Fe	Reserve	2022	Other	4
Santa Fe	Reserve	2022	Other invertebrates	10
Santa Fe	Reserve	2022	Rubble	1
Santa Fe	Reserve	2022	Sand	7
Santa Fe	Reserve	2022	Soft coral	9
Santa Fe	Reserve	2022	Turf algae	19
Santa Fe	Reserve	2023	Hard coral	5
Santa Fe	Reserve	2023	Macroalgae	51
Santa Fe	Reserve	2023	Other	10
Santa Fe	Reserve	2023	Other invertebrates	5
Santa Fe	Reserve	2023	Rubble	3
Santa Fe	Reserve	2023	Sand	7
Santa Fe	Reserve	2023	Soft coral	6
Santa Fe	Reserve	2023	Turf algae	14
Trujillo	Managed Access	2023	Hard coral	4
Trujillo	Managed Access	2023	Macroalgae	66
Trujillo	Managed Access	2023	Other	3
Trujillo	Managed Access	2023	Other invertebrates	4
Trujillo	Managed Access	2023	Sand	6
Trujillo	Managed Access	2023	Soft coral	3
Trujillo	Managed Access	2023	Turf algae	14
Utila	Reserve	2023	Hard coral	9
Utila	Reserve	2023	Macroalgae	17
Utila	Reserve	2023	Other	13
Utila	Reserve	2023	Other invertebrates	7
Utila	Reserve	2023	Rubble	2
Utila	Reserve	2023	Sand	14
Utila	Reserve	2023	Soft coral	6
Utila	Reserve	2023	Turf algae	31

1.2. Mean Benthic Cover Composition per Managed Access Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Normalize and summarize per year and MAA (no management split)
benthic_df_prepped_full <- benthic_summary_combined %>%
  mutate(year = lubridate::year(sample_date)) %>%
  group_by(year, managed_access_areas, category) %>%
  summarise(
    mean_raw_cover = mean(percent_cover, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(year, managed_access_areas) %>%
  mutate(
    total_cover = sum(mean_raw_cover),
    normalized_cover = 100 * mean_raw_cover / total_cover
  ) %>%
  ungroup()

# Step 2: Filter to the most recent year per MAA
latest_years <- benthic_df_prepped_full %>%
  group_by(managed_access_areas) %>%
  summarise(year = max(year), .groups = "drop")

benthic_df_prepped <- benthic_df_prepped_full %>%
  inner_join(latest_years, by = c("managed_access_areas", "year"))

# tep 3: Unique MAAs
maas <- unique(benthic_df_prepped$managed_access_areas)
benthic_composition_stacked <- list()

# Step 4: Define benthic color palette
benthic_colors <- c(
  "Hard coral" = "#C77CFF",
  "Soft coral" = "#F8766D",
  "Macroalgae" = "#66A61E",
  "Turf algae" = "#1B9E77",
  "Other invertebrates" = "#E6AB02",
  "Rubble" = "#A6761D",
  "Sand" = "#FFD700",
  "Bare substrate" = "#666666",
  "Other" = "#999999"
)

# Step 5: Plot per MAA (latest year only)
for (maa in maas) {
  df_subset <- benthic_df_prepped %>%
    filter(managed_access_areas == maa)

  p <- ggplot(df_subset, aes(x = factor(year), y = normalized_cover, fill = category)) +
    geom_bar(stat = "identity", position = "stack", width = 0.7, alpha = 0.8) +
    labs(
      title = paste("Mean Benthic Cover Composition in", maa),
      x = "Year",
      y = "Mean Percent Cover",
      fill = "Benthic Group"
    ) +
    scale_fill_manual(values = benthic_colors) +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "bottom",
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9),
      legend.key.size = unit(0.5, "lines")
    )

  benthic_composition_stacked[[maa]] <- p
}

# Step 6: Display plots
for (maa in maas) {
  print(benthic_composition_stacked[[maa]])
}

# Step 1: Round and select relevant columns (already filtered to latest year)
benthic_table <- benthic_df_prepped %>%
  mutate(normalized_cover = round(normalized_cover, 0)) %>%
  select(managed_access_areas, year, category, normalized_cover)

# Step 2: Format table
table_out <- benthic_table %>%
  kbl(
    caption = "Benthic Cover (%) for Most Recent Year per Managed Access Area",
    col.names = c("Managed Access Area", "Year", "Benthic Category", "Percent Cover"),
    align = "llcr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 3: Add horizontal lines after each Province
break_rows <- benthic_table %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 4: Display final table
table_out

Benthic Cover (%) for Most Recent Year per Managed Access Area

Managed Access Area	Year	Benthic Category	Percent Cover
Iriona and Limon	2021	Bare substrate	1
Iriona and Limon	2021	Hard coral	8
Iriona and Limon	2021	Macroalgae	33
Iriona and Limon	2021	Other	2
Iriona and Limon	2021	Other invertebrates	8
Iriona and Limon	2021	Rubble	1
Iriona and Limon	2021	Sand	35
Iriona and Limon	2021	Soft coral	4
Iriona and Limon	2021	Turf algae	7
Puerto Cortes	2022	Hard coral	3
Puerto Cortes	2022	Macroalgae	24
Puerto Cortes	2022	Other	4
Puerto Cortes	2022	Other invertebrates	3
Puerto Cortes	2022	Sand	8
Puerto Cortes	2022	Soft coral	6
Puerto Cortes	2022	Turf algae	52
Guanaja	2023	Hard coral	9
Guanaja	2023	Macroalgae	39
Guanaja	2023	Other	11
Guanaja	2023	Other invertebrates	1
Guanaja	2023	Rubble	1
Guanaja	2023	Sand	9
Guanaja	2023	Soft coral	10
Guanaja	2023	Turf algae	20
Roatan	2023	Bare substrate	1
Roatan	2023	Hard coral	5
Roatan	2023	Macroalgae	20
Roatan	2023	Other	40
Roatan	2023	Other invertebrates	2
Roatan	2023	Rubble	3
Roatan	2023	Sand	5
Roatan	2023	Soft coral	7
Roatan	2023	Turf algae	16
Santa Fe	2023	Hard coral	5
Santa Fe	2023	Macroalgae	49
Santa Fe	2023	Other	10
Santa Fe	2023	Other invertebrates	8
Santa Fe	2023	Rubble	2
Santa Fe	2023	Sand	7
Santa Fe	2023	Soft coral	5
Santa Fe	2023	Turf algae	13
Trujillo	2023	Hard coral	4
Trujillo	2023	Macroalgae	66
Trujillo	2023	Other	3
Trujillo	2023	Other invertebrates	4
Trujillo	2023	Sand	6
Trujillo	2023	Soft coral	3
Trujillo	2023	Turf algae	14
Utila	2023	Hard coral	9
Utila	2023	Macroalgae	17
Utila	2023	Other	13
Utila	2023	Other invertebrates	7
Utila	2023	Rubble	2
Utila	2023	Sand	14
Utila	2023	Soft coral	6
Utila	2023	Turf algae	31

Mean Benthic Cover by Group

2.1. Mean Benthic Cover by Group by Management Type and Year in Rare Honduras Sites

Interpretation: This bar chart presents the mean percent cover of individual benthic groups averaged across all years and Managed Access Areas. Each colored bar represents the mean contribution of a specific benthic group to total benthic cover under a given management type, with black error bars indicating the standard error. This visualization provides a detailed breakdown of benthic community structure, highlighting which groups are most dominant.

# Normalize and summarize per transect before aggregating
benthic_df_prepped <- benthic_summary_combined %>%
  mutate(year = year(sample_date)) %>%
  group_by(year, management, managed_access_areas, category) %>%
  summarise(
    mean_raw_cover = mean(percent_cover, na.rm = TRUE),
    se_raw_cover = sd(percent_cover, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  ) %>%
  group_by(year, management, managed_access_areas) %>%
  mutate(
    total_cover = sum(mean_raw_cover),
    normalized_cover = 100 * mean_raw_cover / total_cover,
    normalized_se = 100 * se_raw_cover / total_cover
  ) %>%
  ungroup()

# Ensure consistent category order (optional)
benthic_df_prepped$category <- factor(benthic_df_prepped$category)

# List of plots
maas <- unique(benthic_df_prepped$managed_access_areas)
benthic_composition_grouped <- list()

for (maa in maas) {
  df_subset <- benthic_df_prepped %>%
    filter(managed_access_areas == maa)

  benthic_colors <- c(
    "Hard coral" = "#C77CFF",
    "Soft coral" = "#F8766D",
    "Macroalgae" = "#66A61E",
    "Turf algae" = "#1B9E77",
    "Other invertebrates" = "#E6AB02",
    "Rubble" = "#A6761D",
    "Sand" = "#FFD700",
    "Bare substrate" = "#666666",
    "Other" = "#999999"
  )
  
  p <- ggplot(df_subset, aes(x = management, y = normalized_cover, fill = category)) +
    geom_col(position = position_dodge(width = 0.8), width = 0.7, alpha = 0.8) +
    geom_errorbar(
      aes(ymin = normalized_cover - normalized_se, ymax = normalized_cover + normalized_se),
      position = position_dodge(width = 0.8),
      width = 0.2,
      color = "black"
    ) +
    facet_wrap(~ year) +
    labs(
      title = paste("Mean Benthic Cover Composition in", maa),
      x = "Management Type",
      y = "Mean Percent Cover",
      fill = "Benthic Group"
    ) +
    theme_minimal(base_size = 12) +
    scale_fill_manual(values = benthic_colors) +
    theme(
      axis.text.x = element_text(angle = 0),
      legend.position = "bottom",
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9),
      legend.key.size = unit(0.5, "lines")
    )

  benthic_composition_grouped[[maa]] <- p
}

# Print plots
for (maa in maas) {
  print(benthic_composition_grouped[[maa]])
}

2.2. Mean Benthic Cover per Managed Access Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Normalize and summarize across management, for most recent year only

# First, compute full summary by year × MAA × category
benthic_df_prepped_full <- benthic_summary_combined %>%
  mutate(year = lubridate::year(sample_date)) %>%
  group_by(year, managed_access_areas, category) %>%
  summarise(
    mean_raw_cover = mean(percent_cover, na.rm = TRUE),
    se_raw_cover = sd(percent_cover, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  )

# Get most recent year for each MAA
latest_years <- benthic_df_prepped_full %>%
  group_by(managed_access_areas) %>%
  summarise(year = max(year), .groups = "drop")

# Filter and normalize
benthic_df_prepped <- benthic_df_prepped_full %>%
  inner_join(latest_years, by = c("managed_access_areas", "year")) %>%
  group_by(managed_access_areas) %>%
  mutate(
    total_cover = sum(mean_raw_cover),
    normalized_cover = 100 * mean_raw_cover / total_cover,
    normalized_se = 100 * se_raw_cover / total_cover
  ) %>%
  ungroup()

# Optional: Consistent category order
benthic_df_prepped$category <- factor(benthic_df_prepped$category)

# Step 2: Plotting per MAA (only most recent year)
maas <- unique(benthic_df_prepped$managed_access_areas)
benthic_composition_grouped <- list()

benthic_colors <- c(
  "Hard coral" = "#C77CFF",
  "Soft coral" = "#F8766D",
  "Macroalgae" = "#66A61E",
  "Turf algae" = "#1B9E77",
  "Other invertebrates" = "#E6AB02",
  "Rubble" = "#A6761D",
  "Sand" = "#FFD700",
  "Bare substrate" = "#666666",
  "Other" = "#999999"
)

for (maa in maas) {
  df_subset <- benthic_df_prepped %>%
    filter(managed_access_areas == maa)

  p <- ggplot(df_subset, aes(x = category, y = normalized_cover, fill = category)) +
    geom_col(position = position_dodge(width = 0.8), width = 0.7, alpha = 0.8) +
    geom_errorbar(
      aes(ymin = normalized_cover - normalized_se, ymax = normalized_cover + normalized_se),
      position = position_dodge(width = 0.8),
      width = 0.2,
      color = "black"
    ) +
    labs(
      title = paste("Mean Benthic Cover Composition in", maa, "-", unique(df_subset$year )),
      x = "Benthic Group",
      y = "Mean Percent Cover",
      fill = "Benthic Group"
    ) +
    theme_minimal(base_size = 12) +
    scale_fill_manual(values = benthic_colors) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "bottom",
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9),
      legend.key.size = unit(0.5, "lines")
    )

  benthic_composition_grouped[[maa]] <- p
}

# Step 3: Print all plots
for (maa in maas) {
  print(benthic_composition_grouped[[maa]])
}

Mean Hard Coral Cover

3.1. Mean Hard Coral vs. Other Benthic Cover by Management Type and Year in Rare Honduras Sites

Interpretation: This stacked bar chart shows the mean percent cover of hard coral (in purple) compared to all other benthic groups (in gray) across Managed Access and Reserve sites over multiple years and across all Managed Access Areas. Each bar is divided into two segments: the purple portion represents the average hard coral cover, while the gray area represents the combined cover of all other benthic components.

# Step 1: Calculate transect-level cover per original category
benthic_df_detailed <- benthic_summary_combined %>%
  mutate(year = lubridate::year(sample_date)) %>%
  group_by(year, managed_access_areas, management, transect, category) %>%
  summarise(
    transect_cover = mean(percent_cover, na.rm = TRUE),
    .groups = "drop"
  )

# Step 2: Group to "Hard coral" vs "Other" AFTER summarizing per transect
benthic_grouped_se <- benthic_df_detailed %>%
  mutate(category = ifelse(category == "Hard coral", "Hard coral", "Other")) %>%
  group_by(year, managed_access_areas, management, transect, category) %>%
  summarise(transect_cover = sum(transect_cover), .groups = "drop") %>%
  group_by(year, managed_access_areas, management, category) %>%
  summarise(
    mean_cover = mean(transect_cover, na.rm = TRUE),
    se_cover = sd(transect_cover, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  )

# Step 3: Normalize and compute normalized SE
benthic_normalized <- benthic_grouped_se %>%
  group_by(year, managed_access_areas, management) %>%
  mutate(
    total_mean_cover = sum(mean_cover),
    normalized_cover = 100 * mean_cover / total_mean_cover,
    normalized_se = ifelse(category == "Hard coral", 100 * se_cover / total_mean_cover, NA)
  ) %>%
  ungroup() %>%
  mutate(
    category = factor(category, levels = c("Hard coral", "Other"))  # hard coral on bottom
  )

# Step 4: Plot without error bars
maas <- unique(benthic_normalized$managed_access_areas)
benthic_composition_stacked <- list()

benthic_colors <- c(
  "Hard coral" = "#C77CFF",
  "Other" = "#999999"
)

for (maa in maas) {
  df_subset <- benthic_normalized %>%
    filter(managed_access_areas == maa)

  p <- ggplot(df_subset, aes(x = management, y = normalized_cover, fill = category)) +
    geom_bar(
      stat = "identity",
      position = position_stack(reverse = TRUE),  # keeps Hard coral on bottom
      width = 0.7,
      alpha = 0.8
    ) +
    facet_wrap(~ year) +
    labs(
      title = paste("Mean Benthic Cover Composition in", maa),
      x = "Management Type",
      y = "Mean Percent Cover",
      fill = "Benthic Group"
    ) +
    scale_fill_manual(values = benthic_colors) +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "bottom",
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 9),
      legend.key.size = unit(0.5, "lines")
    )

  benthic_composition_stacked[[maa]] <- p
}

# Step 5: Print all plots
for (maa in maas) {
  print(benthic_composition_stacked[[maa]])
}

# Step 1: Create summary table
benthic_table <- benthic_normalized %>%
  group_by(managed_access_areas, management, year, category) %>%
  summarise(normalized_cover = mean(normalized_cover, na.rm = TRUE), .groups = "drop") %>%
  mutate(normalized_cover = round(normalized_cover, 0)) %>%
  select(managed_access_areas, management, year, category, normalized_cover)

# Sep 2: Format table with kableExtra


table_out <- benthic_table %>%
  kbl(
    caption = "Benthic Cover (%) by Managed Access Area, Management Type, and Year (Hard Coral vs. Other)",
    col.names = c("Managed Access Area", "Management Type", "Year", "Benthic Category", "Percent Cover"),
    align = "lllcr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

# Step 3: Add horizontal line after each Province/Management/Year group
break_rows <- benthic_table %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas, management, year) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 4: Show table
table_out

Benthic Cover (%) by Managed Access Area, Management Type, and Year (Hard Coral vs. Other)

Managed Access Area	Management Type	Year	Benthic Category	Percent Cover
Guanaja	Managed Access	2021	Hard coral	6
Guanaja	Managed Access	2021	Other	94
Guanaja	Managed Access	2022	Hard coral	4
Guanaja	Managed Access	2022	Other	96
Guanaja	Reserve	2021	Hard coral	1
Guanaja	Reserve	2021	Other	99
Guanaja	Reserve	2022	Hard coral	4
Guanaja	Reserve	2022	Other	96
Guanaja	Reserve	2023	Hard coral	9
Guanaja	Reserve	2023	Other	91
Iriona and Limon	Managed Access	2021	Hard coral	8
Iriona and Limon	Managed Access	2021	Other	92
Iriona and Limon	Reserve	2021	Hard coral	9
Iriona and Limon	Reserve	2021	Other	91
Puerto Cortes	Reserve	2022	Hard coral	3
Puerto Cortes	Reserve	2022	Other	97
Roatan	Managed Access	2021	Hard coral	15
Roatan	Managed Access	2021	Other	85
Roatan	Managed Access	2023	Hard coral	5
Roatan	Managed Access	2023	Other	95
Roatan	Reserve	2021	Hard coral	16
Roatan	Reserve	2021	Other	84
Santa Fe	Managed Access	2022	Hard coral	9
Santa Fe	Managed Access	2022	Other	91
Santa Fe	Managed Access	2023	Hard coral	8
Santa Fe	Managed Access	2023	Other	92
Santa Fe	Reserve	2021	Hard coral	8
Santa Fe	Reserve	2021	Other	92
Santa Fe	Reserve	2022	Hard coral	13
Santa Fe	Reserve	2022	Other	87
Santa Fe	Reserve	2023	Hard coral	4
Santa Fe	Reserve	2023	Other	96
Trujillo	Managed Access	2023	Hard coral	4
Trujillo	Managed Access	2023	Other	96
Utila	Reserve	2023	Hard coral	9
Utila	Reserve	2023	Other	91

3.2. Mean Hard Coral vs. Other Benthic Cover per Managed Acces Area for the Most Recent Year in Rare Honduras Sites

# Step 1: Calculate transect-level cover per original category
benthic_df_detailed <- benthic_summary_combined %>%
  mutate(year = lubridate::year(sample_date)) %>%
  group_by(year, managed_access_areas, management, transect, category) %>%
  summarise(
    transect_cover = mean(percent_cover, na.rm = TRUE),
    .groups = "drop"
  )

# Step 2: Group to "Hard coral" vs "Other" AFTER summarizing per transect
benthic_grouped_se <- benthic_df_detailed %>%
  mutate(category = ifelse(category == "Hard coral", "Hard coral", "Other")) %>%
  group_by(year, managed_access_areas, management, transect, category) %>%
  summarise(transect_cover = sum(transect_cover), .groups = "drop") %>%
  group_by(year, managed_access_areas, management, category) %>%
  summarise(
    mean_cover = mean(transect_cover, na.rm = TRUE),
    se_cover = sd(transect_cover, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  )

# Step 3: Normalize and compute normalized SE
benthic_normalized <- benthic_grouped_se %>%
  group_by(year, managed_access_areas, management) %>%
  mutate(
    total_mean_cover = sum(mean_cover),
    normalized_cover = 100 * mean_cover / total_mean_cover,
    normalized_se = ifelse(category == "Hard coral", 100 * se_cover / total_mean_cover, NA)
  ) %>%
  ungroup() %>%
  mutate(
    category = factor(category, levels = c("Hard coral", "Other"))
  )

# Step 4: Get latest year per MAA
latest_years <- benthic_normalized %>%
  group_by(managed_access_areas) %>%
  summarise(latest_year = max(year), .groups = "drop")

# Step 5: Filter to latest year and drop management
benthic_latest <- benthic_normalized %>%
  inner_join(latest_years, by = "managed_access_areas") %>%
  filter(year == latest_year) %>%
  group_by(managed_access_areas, year, category) %>%
  summarise(
    normalized_cover = mean(normalized_cover, na.rm = TRUE),
    .groups = "drop"
  )

# Step 6: Plot (1 plot per MAA with year on x-axis and removed from title)
benthic_colors <- c(
  "Hard coral" = "#C77CFF",
  "Other" = "#999999"
)

maas <- unique(benthic_latest$managed_access_areas)
benthic_composition_latest <- list()

for (maa in maas) {
  df_subset <- benthic_latest %>%
    filter(managed_access_areas == maa)

  p <- ggplot(df_subset, aes(x = factor(year), y = normalized_cover, fill = category)) +
    geom_bar(stat = "identity", position = position_stack(reverse = TRUE), width = 0.6, alpha = 0.8) +
    labs(
      title = paste("Mean Benthic Cover in", maa),
      x = "Year",
      y = "Mean Percent Cover",
      fill = "Benthic Group"
    ) +
    scale_fill_manual(values = benthic_colors) +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "bottom"
    )

  benthic_composition_latest[[maa]] <- p
}

# Step 7: Print plots
for (maa in maas) {
  print(benthic_composition_latest[[maa]])
}

# Step 1: Create summary table
benthic_table_latest <- benthic_latest %>%
  mutate(normalized_cover = round(normalized_cover, 0)) %>%
  select(managed_access_areas, year, category, normalized_cover)

# Step 2: Format table with horizontal lines by province
break_rows <- benthic_table_latest %>%
  mutate(row_number = row_number()) %>%
  group_by(managed_access_areas) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

table_out <- benthic_table_latest %>%
  kbl(
    caption = "Benthic Cover (%) for Most Recent Year per Managed Access Area",
    col.names = c("Managed Access Area", "Year", "Benthic Category", "Percent Cover"),
    align = "llcr",
    booktabs = TRUE
  ) %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))

for (row in break_rows) {
  table_out <- table_out %>%
    row_spec(row, extra_css = "border-bottom: 2px solid #000;")
}

# Step 3: Show table
table_out

Benthic Cover (%) for Most Recent Year per Managed Access Area

Managed Access Area	Year	Benthic Category	Percent Cover
Guanaja	2023	Hard coral	9
Guanaja	2023	Other	91
Iriona and Limon	2021	Hard coral	8
Iriona and Limon	2021	Other	92
Puerto Cortes	2022	Hard coral	3
Puerto Cortes	2022	Other	97
Roatan	2023	Hard coral	5
Roatan	2023	Other	95
Santa Fe	2023	Hard coral	6
Santa Fe	2023	Other	94
Trujillo	2023	Hard coral	4
Trujillo	2023	Other	96
Utila	2023	Hard coral	9
Utila	2023	Other	91