Data Wrangling

We created these aggregated hard-coral datasets by integrating three data sources: OP41 data, historical Fish Forever 2.0 data, and MERMAID data. Before combining them, we carried out a duplicate review across sources and removed overlapping records so the same surveys would not be counted twice. We then estimated hard coral cover (%) at the sample level using a source-specific approach that preserved each dataset’s native structure, harmonized clearly equivalent ABM names across datasets, and excluded records with missing ABM assignments. OP41 Indonesia was left out of the integrated hard-coral analysis because that dataset did not include the metadata needed to match listed sites to ABMs. We also excluded observations where management regime could not be confidently identified, because some MERMAID records had management information that was not entered in a usable or standardized way. Finally, we aggregated the valid sample-level estimates by country, ma_name (ABM), management_name, and sample_year, while excluding samples that failed basic quality-control checks.

# ============================================================
# Build harmonized hard-coral datasets from raw benthic sources
#
# Source families integrated:
# - OP41 data (Philippines only in final hard-coral integration)
# - Historical Fish Forever 2.0 (FF2.0) data
# - MERMAID data (PIT and PQT)
#
# Main steps:
# 1. Read raw source datasets and lookup tables
# 2. Remove known duplicate records across sources
# 3. Add / standardize ma_name and management_name where needed
# 4. Harmonize a few clearly equivalent ABM names
# 5. Estimate hard coral cover (%) at the sample level
# 6. Aggregate valid samples by country + ma_name + management_name + year
#
# Final outputs:
# - hard_coral_by_ma_mgmt_year
# - hard_coral_samples_export
# - hard_coral_samples_excluded_qc
#
# Important:
# - OP41 Indonesia is read and partially cleaned for duplicate control,
#   but it is not included in the final hard-coral integration because
#   that dataset does not include metadata to match site codes to ABMs.
# ============================================================


# ------------------------------------------------------------
# 1) Read raw source datasets
# ------------------------------------------------------------

mermaid_pqt <- read_csv(
  here("data", "sql_mv_benthic_pqt_sum_2026-03-27_103440.csv"),
  show_col_types = FALSE
)

mermaid_pit <- read_csv(
  here("data", "sql_mv_benthic_pit_sum_2026-03-27_103045.csv"),
  show_col_types = FALSE
)

ff2_habitat <- read_csv(
  here("data", "ff2_habitat.csv"),
  show_col_types = FALSE
)

op41_habitat_phl <- read_csv(
  here("data", "op41_habitat_phl.csv"),
  show_col_types = FALSE
)

# OP41 Indonesia is not UTF-8, so read it as Latin1
op41_habitat_ind <- read_csv(
  here("data", "op41_habitat_ind.csv"),
  locale = locale(encoding = "Latin1"),
  show_col_types = FALSE
)

# ------------------------------------------------------------
# 2) Read lookup tables used to assign ABMs (ma_name)
# ------------------------------------------------------------

op41_site_codes <- read_excel(
  here("data", "op41_site_codes.xlsx")
)

hon_sites <- read_excel(
  here("data", "hon_ff2.0_habitat_sites_processed.xlsx")
)

phi_sites <- read_excel(
  here("data", "phi_ff2.0_habitat_sites_processed.xlsx")
)

ind_sites <- read_excel(
  here("data", "ind_ff2.0_habitat_sites_processed.xlsx")
)

# ------------------------------------------------------------
# 3) Remove known duplicate records across sources
# ------------------------------------------------------------

# 3a. OP41 Indonesia: remove KOLONOB 2017 because those records
#     are already represented in FF2.0
op41_habitat_ind_original_n <- nrow(op41_habitat_ind)

op41_habitat_ind <- op41_habitat_ind %>%
  filter(!(site_code == "KOLONOB" & Year == 2017))

op41_removed_n <- op41_habitat_ind_original_n - nrow(op41_habitat_ind)

# 3b. MERMAID PIT: remove Indonesia 2021 block because it is
#     already represented in FF2.0
mermaid_pit_original_n <- nrow(mermaid_pit)

mermaid_pit <- mermaid_pit %>%
  filter(!(
    country_name == "Indonesia" &
      year == 2021 &
      project_name == "SE Sulawesi Monitoring Biofisik 2021"
  ))

pit_removed_n <- mermaid_pit_original_n - nrow(mermaid_pit)

# 3c. FF2.0: remove selected Honduras site-year duplicates where
#     MERMAID is treated as the preferred source
honduras_dup_sites_2021 <- c("6", "7", "8", "10", "11", "17", "18", "19", "20", "21", "22")
honduras_dup_sites_2023 <- c("UOCA1", "UOIA2", "UWIA2")

ff2_habitat_original_n <- nrow(ff2_habitat)

ff2_habitat <- ff2_habitat %>%
  filter(!(
    country == "Honduras" &
      methodology == "photo transect" &
      (
        (year == 2021 & location_name %in% honduras_dup_sites_2021) |
        (year == 2023 & location_name %in% honduras_dup_sites_2023)
      )
  ))

ff2_removed_n <- ff2_habitat_original_n - nrow(ff2_habitat)

# ------------------------------------------------------------
# 4) Add / standardize ma_name and management_name
# ------------------------------------------------------------

# 4a. Build FF2.0 lookup to assign ma_name from processed site files
ff2_sites_lookup <- bind_rows(hon_sites, phi_sites, ind_sites) %>%
  transmute(
    country = str_trim(country),
    join_name = name %>%
      str_trim() %>%
      str_to_lower() %>%
      str_replace_all(" ", "_"),
    ma_name = site_name %>% str_trim()
  ) %>%
  distinct()

# 4b. Add ma_name to FF2.0
ff2_habitat <- ff2_habitat %>%
  select(-any_of(c("ma_name", "ma_name.x", "ma_name.y"))) %>%
  mutate(
    country = str_trim(country),
    join_name = location_name %>%
      str_trim() %>%
      str_to_lower() %>%
      str_replace_all(" ", "_")
  ) %>%
  left_join(ff2_sites_lookup, by = c("country", "join_name")) %>%
  select(-join_name) %>%
  mutate(
    # where lookup is missing, fall back to sitename if available
    ma_name = coalesce(ma_name, sitename)
  )

# 4c. Add ma_name to OP41 Philippines
op41_phl_lookup <- op41_site_codes %>%
  transmute(
    site_code = `Campaign Site Abbreviated Name` %>% str_trim(),
    ma_name = `Campaign Site Full Name` %>% str_trim()
  ) %>%
  distinct()

op41_habitat_phl <- op41_habitat_phl %>%
  select(-any_of(c("ma_name", "ma_name.x", "ma_name.y"))) %>%
  mutate(site_code = str_trim(site_code)) %>%
  left_join(op41_phl_lookup, by = "site_code")

# 4d. Standardize management_name in OP41 Philippines
op41_habitat_phl <- op41_habitat_phl %>%
  rename(management_name = inside_or_outside_NTZ) %>%
  mutate(
    management_name = case_when(
      str_to_lower(management_name) == "inside" ~ "Reserve",
      str_to_lower(management_name) %in% c("outside", "turf") ~ "Managed Access",
      TRUE ~ management_name
    )
  )

# 4e. Standardize management_name in OP41 Indonesia
#     (kept for reference / duplicate control, but not integrated)
op41_habitat_ind <- op41_habitat_ind %>%
  rename(management_name = inside_or_outside_NTZ) %>%
  mutate(
    management_name = case_when(
      str_to_lower(management_name) == "inside" ~ "Reserve",
      str_to_lower(management_name) %in% c("outside", "turf") ~ "Managed Access",
      TRUE ~ management_name
    )
  )

# 4f. Standardize management_name in FF2.0
ff2_habitat <- ff2_habitat %>%
  rename(management_name = location_status)

# ------------------------------------------------------------
# 5) Harmonize a few clearly equivalent ABM names
# ------------------------------------------------------------

harmonize_ma_name <- function(x) {
  x <- str_squish(x)

  case_when(
    x == "Puerto Cortés" ~ "Puerto Cortes",
    x == "Iriona and Limón" ~ "Iriona and Limon",
    x == "San Carlos City" ~ "San Carlos",
    x == "Sta. Monica" ~ "Santa Monica",
    TRUE ~ x
  )
}

op41_habitat_phl <- op41_habitat_phl %>%
  mutate(ma_name = harmonize_ma_name(ma_name))

ff2_habitat <- ff2_habitat %>%
  mutate(ma_name = harmonize_ma_name(ma_name))

mermaid_pit <- mermaid_pit %>%
  mutate(ma_name = harmonize_ma_name(ma_name))

mermaid_pqt <- mermaid_pqt %>%
  mutate(ma_name = harmonize_ma_name(ma_name))

# ------------------------------------------------------------
# 6) Build sample-level hard coral cover
# ------------------------------------------------------------

# 6a. OP41 Philippines
# Sample = one transect_id on one survey date
# Hard coral cover = percent_live_coral_cover
op41_hard_coral_samples <- op41_habitat_phl %>%
  filter(!is.na(ma_name), ma_name != "") %>%
  mutate(
    source = "op41_habitat_phl",
    survey_date = suppressWarnings(mdy(data_collection_date)),
    sample_year = coalesce(year(survey_date), as.integer(Year)),
    sample_unit_id = as.character(transect_id),
    hard_coral_cover_pct = percent_live_coral_cover,
    sample_weight = 100,
    total_cover_check =
      percent_live_coral_cover +
      percent_soft_coral_cover +
      percent_algae_cover +
      percent_other_biota_cover +
      percent_available_substrate_cover +
      percent_mobile_substrate_cover +
      percent_BLK,
    use_in_summary = between(total_cover_check, 90, 110)
  ) %>%
  select(
    source, country, ma_name, management_name,
    survey_date, sample_year, sample_unit_id,
    hard_coral_cover_pct, sample_weight,
    total_cover_check, use_in_summary
  )

# 6b. FF2.0
# Sample = one transect_no at one site on one survey date
# Hard coral cover = sum of percentage where category == "Hard coral"
ff2_hard_coral_samples <- ff2_habitat %>%
  distinct() %>%
  filter(!is.na(ma_name), ma_name != "") %>%
  mutate(
    survey_date = suppressWarnings(ymd(surveydate)),
    sample_year = coalesce(year(survey_date), as.integer(year)),
    site_key = case_when(
      !is.na(lat) & !is.na(lon) ~ paste0(round(lat, 5), "__", round(lon, 5)),
      !is.na(location_name) & location_name != "" ~ str_to_lower(str_squish(location_name)),
      !is.na(sitename) & sitename != "" ~ str_to_lower(str_squish(sitename)),
      TRUE ~ NA_character_
    )
  ) %>%
  group_by(
    country, ma_name, management_name,
    survey_date, sample_year, site_key, transect_no
  ) %>%
  summarise(
    source = "ff2_habitat",
    hard_coral_cover_pct = sum(percentage[category == "Hard coral"], na.rm = TRUE),
    total_cover_check = sum(percentage, na.rm = TRUE),
    sample_weight = 100,
    .groups = "drop"
  ) %>%
  mutate(
    sample_unit_id = paste(
      country, ma_name, management_name,
      coalesce(as.character(survey_date), as.character(sample_year)),
      site_key, transect_no,
      sep = "__"
    ),
    use_in_summary = between(total_cover_check, 90, 110)
  ) %>%
  select(
    source, country, ma_name, management_name,
    survey_date, sample_year, sample_unit_id,
    hard_coral_cover_pct, sample_weight,
    total_cover_check, use_in_summary
  )

# 6c. MERMAID PIT
# Sample = one site_id on one survey date
# Hard coral cover = hard coral point_count / total point_count * 100
pit_hard_coral_samples <- mermaid_pit %>%
  filter(!is.na(ma_name), ma_name != "") %>%
  mutate(
    survey_date = suppressWarnings(ymd(date)),
    sample_year = coalesce(year(survey_date), as.integer(year))
  ) %>%
  group_by(
    country = country_name, ma_name, management_name,
    survey_date, sample_year, site_id
  ) %>%
  summarise(
    source = "mermaid_pit",
    hard_coral_n = sum(point_count[benthic_category == "Hard coral"], na.rm = TRUE),
    total_n = sum(point_count, na.rm = TRUE),
    hard_coral_cover_pct = if_else(total_n > 0, 100 * hard_coral_n / total_n, NA_real_),
    sample_weight = total_n,
    total_cover_check = total_n,
    .groups = "drop"
  ) %>%
  mutate(
    sample_unit_id = paste(site_id, survey_date, sep = "__"),
    use_in_summary = total_cover_check > 0 & !is.na(hard_coral_cover_pct)
  ) %>%
  select(
    source, country, ma_name, management_name,
    survey_date, sample_year, sample_unit_id,
    hard_coral_cover_pct, sample_weight,
    total_cover_check, use_in_summary
  )

# 6d. MERMAID PQT
# Sample = one site_id in one year
# Hard coral cover = hard coral observation_count / total observation_count * 100
# Note: this source has no exact survey date
pqt_hard_coral_samples <- mermaid_pqt %>%
  filter(!is.na(ma_name), ma_name != "") %>%
  group_by(
    country = country_name, ma_name, management_name,
    year, site_id
  ) %>%
  summarise(
    source = "mermaid_pqt",
    hard_coral_n = sum(observation_count[benthic_category == "Hard coral"], na.rm = TRUE),
    total_n = sum(observation_count, na.rm = TRUE),
    hard_coral_cover_pct = if_else(total_n > 0, 100 * hard_coral_n / total_n, NA_real_),
    sample_weight = total_n,
    total_cover_check = total_n,
    .groups = "drop"
  ) %>%
  mutate(
    survey_date = as.Date(NA),
    sample_year = as.integer(year),
    sample_unit_id = paste(site_id, sample_year, sep = "__"),
    use_in_summary = total_cover_check > 0 & !is.na(hard_coral_cover_pct)
  ) %>%
  select(
    source, country, ma_name, management_name,
    survey_date, sample_year, sample_unit_id,
    hard_coral_cover_pct, sample_weight,
    total_cover_check, use_in_summary
  )

# ------------------------------------------------------------
# 7) Combine all sample-level hard coral datasets
# ------------------------------------------------------------

hard_coral_samples <- bind_rows(
  op41_hard_coral_samples,
  ff2_hard_coral_samples,
  pit_hard_coral_samples,
  pqt_hard_coral_samples
)

# ------------------------------------------------------------
# 8) Create final output datasets
# ------------------------------------------------------------

# 8a. Main yearly summary
hard_coral_by_ma_mgmt_year <- hard_coral_samples %>%
  filter(use_in_summary) %>%
  group_by(country, ma_name, management_name, sample_year) %>%
  summarise(
    n_samples = n(),
    n_sources = n_distinct(source),
    sources = paste(sort(unique(source)), collapse = "; "),
    mean_hard_coral_cover_pct = mean(hard_coral_cover_pct, na.rm = TRUE),
    weighted_hard_coral_cover_pct = weighted.mean(
      hard_coral_cover_pct,
      w = sample_weight,
      na.rm = TRUE
    ),
    .groups = "drop"
  ) %>%
  arrange(country, ma_name, management_name, sample_year)

# 8b. Sample-level export used in the summary
hard_coral_samples_export <- hard_coral_samples %>%
  filter(use_in_summary) %>%
  select(
    country,
    ma_name,
    management_name,
    survey_date,
    sample_year,
    hard_coral_cover_pct,
    source
  ) %>%
  arrange(country, ma_name, management_name, sample_year, survey_date, source)

# 8c. Sample-level QC table of excluded rows
hard_coral_samples_excluded_qc <- hard_coral_samples %>%
  filter(!use_in_summary) %>%
  arrange(country, ma_name, management_name, sample_year, survey_date, source)

# ------------------------------------------------------------
# 9) Messages and quick inspection
# ------------------------------------------------------------

message("Rows removed from op41_habitat_ind (KOLONOB 2017): ", op41_removed_n)
message("Rows removed from MERMAID PIT (Indonesia 2021 duplicates): ", pit_removed_n)
message("Rows removed from ff2_habitat (Honduras 2021 + 2023 duplicates): ", ff2_removed_n)

message("Rows in hard_coral_samples: ", nrow(hard_coral_samples))
message("Rows used in yearly summary: ", sum(hard_coral_samples$use_in_summary, na.rm = TRUE))
message("Rows excluded by QC: ", sum(!hard_coral_samples$use_in_summary, na.rm = TRUE))

message(
  "Rows dropped because ma_name was missing in MERMAID PQT: ",
  sum(is.na(mermaid_pqt$ma_name) | mermaid_pqt$ma_name == "", na.rm = TRUE)
)

message(
  "Unmatched ff2_habitat rows after adding ma_name: ",
  ff2_habitat %>% filter(is.na(ma_name)) %>% nrow()
)

message(
  "Unmatched op41_habitat_phl rows after adding ma_name: ",
  op41_habitat_phl %>% filter(is.na(ma_name)) %>% nrow()
)


# Remove non standard management
hard_coral_samples_export <- hard_coral_samples_export %>%
  filter(management_name %in% c("Reserve", "Managed Access"))


# ------------------------------------------------------------
# 10) Main outputs
# ------------------------------------------------------------
# hard_coral_by_ma_mgmt_year # Main summary table. It was built by first estimating hard coral cover % at the sample level for each source dataset, then combining all valid samples from the 4 sources, and finally summarizing them by year
# hard_coral_samples_export # Sample-level dataset used for hard_coral_by_ma_mgmt_year. 
# hard_coral_samples_excluded_qc # Sample-level QC table of rows that were excluded from the summary. It is useful for checking problematic samples.

# ------------------------------------------------------------
# 11) Optional write-out
# ------------------------------------------------------------

# write_csv(hard_coral_by_ma_mgmt_year, here("data", "hard_coral_by_ma_mgmt_year.csv"))
# write_csv(hard_coral_samples_export, here("data", "hard_coral_samples_export.csv"))
# write_csv(hard_coral_samples_excluded_qc, here("data", "hard_coral_samples_excluded_qc.csv"))

Exploratory Analysis

We used hard coral cover because it is the most widely used and comparable indicators of reef condition across monitoring programs. Hard coral cover reflects the abundance of the main reef-building corals and is commonly used as a standard indicator of coral reef health, whereas broader alternatives such as total coral cover are less consistently reported across datasets and can be less comparable across programs.


# ============================================================
# Read analysis-ready sample-level dataset
# ============================================================

hc <- read_csv(
  here("data", "hard_coral_samples_export.csv"),
  show_col_types = FALSE
) %>%
  mutate(
    survey_date = as.Date(survey_date),
    sample_year = as.integer(sample_year),
    management_name = factor(
      management_name,
      levels = c("Managed Access", "Reserve")
    ),
    management_name_safe = recode(
      as.character(management_name),
      "Managed Access" = "Managed_Access",
      "Reserve" = "Reserve"
    ),
    source = factor(source)
  )

# ============================================================
# 1) Table: number of samples per management regime by country
#    ordered by total number of samples
# ============================================================

samples_per_management_country <- hc %>%
  count(country, management_name, name = "n_samples") %>%
  pivot_wider(
    names_from = management_name,
    values_from = n_samples,
    values_fill = 0
  ) %>%
  mutate(total_samples = `Managed Access` + Reserve) %>%
  arrange(desc(total_samples)) %>%
  select(country, `Managed Access`, Reserve, total_samples)

samples_per_management_country
#> # A tibble: 4 × 4
#>   country     `Managed Access` Reserve total_samples
#>   <chr>                  <int>   <int>         <int>
#> 1 Philippines              518     578          1096
#> 2 Indonesia                549     333           882
#> 3 Honduras                  31      93           124
#> 4 Mozambique                 2      12            14

samples_per_management_country %>%
  kbl(
    caption = "Number of hard coral samples by country and management regime",
    align = c("l", "r", "r", "r")
  ) %>%
  kable_styling(full_width = FALSE, position = "left")
Number of hard coral samples by country and management regime
country Managed Access Reserve total_samples
Philippines 518 578 1096
Indonesia 549 333 882
Honduras 31 93 124
Mozambique 2 12 14 

# ============================================================
# 2) Summary statistics for the whole dataset
#    rounded to 1 decimal
# ============================================================

overall_summary <- hc %>%
  summarise(
    n = n(),
    mean = mean(hard_coral_cover_pct, na.rm = TRUE),
    sd = sd(hard_coral_cover_pct, na.rm = TRUE),
    median = median(hard_coral_cover_pct, na.rm = TRUE),
    q1 = quantile(hard_coral_cover_pct, 0.25, na.rm = TRUE),
    q3 = quantile(hard_coral_cover_pct, 0.75, na.rm = TRUE),
    min = min(hard_coral_cover_pct, na.rm = TRUE),
    max = max(hard_coral_cover_pct, na.rm = TRUE)
  ) %>%
  mutate(across(where(is.double), ~ round(.x, 1)))

overall_summary
#> # A tibble: 1 × 8
#>       n  mean    sd median    q1    q3   min   max
#>   <int> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1  2116  34.4  21.3   32.5    17    49     0   100

overall_summary %>%
  kbl(
    caption = "Summary statistics for hard coral cover across the full dataset",
    align = "r"
  ) %>%
  kable_styling(full_width = FALSE, position = "left")
Summary statistics for hard coral cover across the full dataset
n mean sd median q1 q3 min max
2116 34.4 21.3 32.5 17 49 0 100 

# ============================================================
# Build ABM-year summary used by plots below
# ============================================================

abm_year_summary <- hc %>%
  group_by(country, ma_name, management_name, management_name_safe, sample_year) %>%
  summarise(
    n_samples = n(),
    mean_hard_coral_cover_pct = mean(hard_coral_cover_pct, na.rm = TRUE),
    median_hard_coral_cover_pct = median(hard_coral_cover_pct, na.rm = TRUE),
    .groups = "drop"
  )

# ============================================================
# 3) Reserve vs Managed Access within the same ABM-year
#    countries stacked vertically
# ============================================================

reserve_vs_managed_pairs <- abm_year_summary %>%
  filter(country != "Mozambique") %>%
  select(
    country, ma_name, sample_year, management_name_safe,
    n_samples, mean_hard_coral_cover_pct
  ) %>%
  pivot_wider(
    names_from = management_name_safe,
    values_from = c(n_samples, mean_hard_coral_cover_pct),
    names_sep = "__"
  ) %>%
  filter(
    !is.na(mean_hard_coral_cover_pct__Reserve),
    !is.na(mean_hard_coral_cover_pct__Managed_Access)
  ) %>%
  mutate(
    total_samples_pair = n_samples__Reserve + n_samples__Managed_Access,
    diff_reserve_minus_managed =
      mean_hard_coral_cover_pct__Reserve - mean_hard_coral_cover_pct__Managed_Access
  ) %>%
  arrange(country, ma_name, sample_year)

p_reserve_vs_managed_vertical <- ggplot(
  reserve_vs_managed_pairs,
  aes(
    x = mean_hard_coral_cover_pct__Managed_Access,
    y = mean_hard_coral_cover_pct__Reserve
  )
) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  geom_point(aes(size = total_samples_pair), alpha = 0.8) +
  facet_wrap(~ country, scales = "free", ncol = 1) +
  labs(
    title = "Reserve vs Managed Access within the same ABM-year",
    x = "Managed Access mean hard coral cover (%)",
    y = "Reserve mean hard coral cover (%)",
    size = "Total samples\nin pair"
  ) +
  theme_minimal(base_size = 14)

p_reserve_vs_managed_vertical


# ============================================================
# 4) One combined ABM-level trajectory plot
#    - all ABMs in the same plot
#    - countries stacked vertically
#    - thin ABM lines
#    - thicker overall trend lines by management regime
#
# NOTE:
# The trend lines are descriptive smooths across ABM-year means,
# not inferential trend estimates.
# ============================================================

abm_year_summary_nomoz <- abm_year_summary %>%
  filter(country %in% c("Honduras", "Indonesia", "Philippines"))

all_years <- seq(
  min(abm_year_summary_nomoz$sample_year, na.rm = TRUE),
  max(abm_year_summary_nomoz$sample_year, na.rm = TRUE),
  by = 1
)

p_abm_trajectories_combined <- ggplot(
  abm_year_summary_nomoz,
  aes(
    x = sample_year,
    y = mean_hard_coral_cover_pct,
    color = management_name
  )
) +
  # thin ABM-specific trajectories
  geom_line(
    aes(group = interaction(ma_name, management_name)),
    alpha = 0.18,
    linewidth = 0.4
  ) +
  geom_point(
    aes(size = n_samples),
    alpha = 0.30
  ) +
  # overall descriptive trend by management regime within country
  geom_smooth(
    aes(group = management_name),
    method = "loess",
    se = FALSE,
    linewidth = 1
  ) +
  facet_wrap(~ country, scales = "free_y", ncol = 1) +
  scale_x_continuous(
    breaks = all_years,
    minor_breaks = NULL
  ) +
  labs(
    title = "ABM-level hard coral cover trajectories",
    subtitle = "Thin lines show ABM-specific trajectories; thick lines show descriptive trends by management regime",
    x = "Year",
    y = "Mean hard coral cover (%)",
    color = "Management",
    size = "n samples"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

p_abm_trajectories_combined

The integrated dataset contains 2,116 sample-level hard coral observations. Sample coverage is highly uneven across countries: the Philippines (1,096 samples) and Indonesia (882 samples) account for most of the data, while Honduras (124 samples) is much smaller and Mozambique (14 samples) is too limited for meaningful interpretation. Overall, hard coral cover is moderately variable, with a mean of 34.4%, median of 32.5%, and an interquartile range from 17% to 49%, indicating substantial heterogeneity across sites, years, and management areas.

The paired comparison of Reserve vs Managed Access within the same ABM-year suggests that the relationship between management regime and hard coral cover is not uniform across countries. In the Philippines, most paired observations cluster close to the 1:1 line but tend to fall slightly above it, suggesting a modest tendency for reserves to have higher hard coral cover than managed access areas in the same ABM-year. In Indonesia, the same positive relationship is visible, but the points are more widely dispersed, indicating greater heterogeneity and a weaker, less consistent reserve effect. In Honduras, the number of ABM-year pairs is small, so patterns should be interpreted cautiously.

The ABM-level trajectory plots reinforce this interpretation. In the Philippines, hard coral cover appears relatively high in the earlier years and lower in the later years, but this pattern should be interpreted carefully because it may reflect changes in source coverage and monitored ABMs over time rather than a true ecological decline. In Indonesia, coral cover is generally moderate to high, with both reserves and managed access areas showing broad overlap and no strong consistent separation through time. In Honduras, values are much lower overall, and the small number of observations limits inference.

Taken together, these exploratory results suggest three main points:

  1. first, hard coral cover varies strongly across countries and ABMs;

  2. second, the clearest indication of a reserve advantage appears in the Philippines, though the effect is modest;

  3. and third, temporal patterns should be interpreted cautiously because the dataset is unbalanced across years, countries, and source types.

These results are therefore best treated as descriptive evidence that could help identify priorities for more focused follow-up analysis. A useful next step would be to combine hard coral cover with other ecological, fisheries, and socioeconomic indicators (such as fish biomass, catch data, and household or community-level data) to explore relationships across datasets and see whether bringing these different dimensions together helps reveal a clearer picture of what is happening at the ABM level.