FF Brazil: Complementary Oyster Reef Monitoring Analyses
Mariano Viz
2026-06-19
This doc provides complementary analyses/plots requested for the ecological monitoring presentation of oyster reefs. It focuses on four practical questions:
- How does mean oyster density vary by reef and year?
- Is there enough sampling coverage to discuss seasonal patterns between rainy and dry periods?
- How do live/dead oyster proportions vary by reef, year, and season?
- What can the associated biodiversity data suggest about reef condition, especially in relation to oyster density and size structure?
Density interpretation note. Density is calculated
from live + dead oysters. I set the default parameter
quadrat_area_m2 <- 1, which means results are
interpreted as oysters per quadrat. We could show
results as oysters per m² ff the real quadrat size is
known. I just have to update that value in the analysis code. For
example, if quadrats are 50 cm × 50 cm, I set
quadrat_area_m2 <- 0.25 to express density as
oysters per m².
Biodiversity note. The 2021–2024 data includes
associated taxa counts in taxonomic_group_counts. The 2025
data does not have the same associated biodiversity count columns.
Therefore, density and live/dead analyses use 2021–2025 data, while
biodiversity/richness analyses are based on 2021–2024 data only.
# Run analysis code:
# FF Brazil oyster reefs: complementary analyses for final ecological monitoring presentation
# Density, live/dead proportions, seasonality, associated biodiversity, and condition relationships
#
# Designed to live in the same project as the previous Rmds.
# Expected input paths:
# data/raw/BRA_oyster_monitoring_database_2021_2024.xlsx
# data/raw/BRA_oyster_monitoring_database_2025.xlsx
#
# Main outputs:
# outputs/oyster_complementary_analysis/plots/*.png
# outputs/oyster_complementary_analysis/tables/*.csv
# ------------------------------------------------------------------------------
# 0. Packages and setup
# ------------------------------------------------------------------------------
library(tidyverse)
library(readxl)
library(here)
library(janitor)
library(lubridate)
library(scales)
library(kableExtra)
# ---- User-adjustable parameters ----
# IMPORTANT: If the quadrat area is known, update this value.
# - If quadrats are 50 cm x 50 cm, use 0.25.
# - If you leave it as 1, density is interpreted as oysters per quadrat.
quadrat_area_m2 <- 1
# Minimum quadrats used to flag whether reef-year-season comparisons are robust enough
# for presentation interpretation. Adjust if the monitoring protocol uses a different rule.
min_quadrats_main <- 10
min_quadrats_seasonal <- 20
# Output folders
output_dir <- here("outputs", "oyster_complementary_analysis")
plot_dir <- file.path(output_dir, "plots")
table_dir <- file.path(output_dir, "tables")
dir.create(plot_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(table_dir, recursive = TRUE, showWarnings = FALSE)
old_file <- here("data", "raw", "BRA_oyster_monitoring_database_2021_2024.xlsx")
new_file <- here("data", "raw", "BRA_oyster_monitoring_database_2025.xlsx")
# Plot theme used across all presentation charts
theme_oyster <- function(base_size = 12) {
theme_minimal(base_size = base_size) +
theme(
plot.title = element_text(face = "bold", size = base_size + 3),
plot.subtitle = element_text(size = base_size, color = "grey30"),
plot.caption = element_text(size = base_size - 2, color = "grey40", hjust = 0),
panel.grid.minor = element_blank(),
strip.text = element_text(face = "bold"),
legend.position = "bottom"
)
}
save_plot <- function(plot, filename, width = 12, height = 8, dpi = 320) {
ggsave(
filename = file.path(plot_dir, filename),
plot = plot,
width = width,
height = height,
dpi = dpi,
bg = "white"
)
}
parse_numeric_safely <- function(x) {
x_chr <- str_squish(as.character(x))
x_chr <- na_if(x_chr, "")
x_chr <- na_if(x_chr, "NA")
suppressWarnings(as.numeric(x_chr))
}
to_sampling_date <- function(x) {
if (inherits(x, "Date")) return(x)
if (inherits(x, "POSIXt")) return(as.Date(x))
x_chr <- str_squish(as.character(x))
x_chr <- na_if(x_chr, "")
x_chr <- na_if(x_chr, "NA")
x_num <- suppressWarnings(as.numeric(x_chr))
out <- suppressWarnings(as.Date(x_chr))
out[!is.na(x_num)] <- as.Date(x_num[!is.na(x_num)], origin = "1899-12-30")
out
}
normalize_reef <- function(x) {
x_clean <- x %>%
as.character() %>%
str_squish() %>%
na_if("")
case_when(
x_clean == "Terra amarela" ~ "Terra Amarela",
TRUE ~ x_clean
)
}
normalize_season <- function(x) {
x_clean <- x %>% as.character() %>% str_squish() %>% str_to_lower()
case_when(
x_clean %in% c("rainy", "chuvoso", "chuva", "wet") ~ "Rainy",
x_clean %in% c("dry", "seco") ~ "Dry",
TRUE ~ NA_character_
)
}
wilson_ci_low <- function(success, total, z = qnorm(0.975)) {
n <- ifelse(total <= 0 | is.na(total), NA_real_, total)
p <- success / n
denom <- 1 + z^2 / n
out <- ((p + z^2 / (2 * n)) - z * sqrt((p * (1 - p) + z^2 / (4 * n)) / n)) / denom
pmax(0, out)
}
wilson_ci_high <- function(success, total, z = qnorm(0.975)) {
n <- ifelse(total <= 0 | is.na(total), NA_real_, total)
p <- success / n
denom <- 1 + z^2 / n
out <- ((p + z^2 / (2 * n)) + z * sqrt((p * (1 - p) + z^2 / (4 * n)) / n)) / denom
pmin(1, out)
}
safe_cor <- function(x, y, method = "spearman") {
ok <- complete.cases(x, y)
if (sum(ok) < 3) return(NA_real_)
cor(x[ok], y[ok], method = method)
}
shannon_index <- function(counts) {
counts <- counts[!is.na(counts) & counts > 0]
if (length(counts) == 0 || sum(counts) == 0) return(NA_real_)
p <- counts / sum(counts)
-sum(p * log(p))
}
density_unit_label <- if (quadrat_area_m2 == 1) {
"oysters per quadrat"
} else {
"oysters per m²"
}
# ------------------------------------------------------------------------------
# 1. Read data using the same approach as the existing Rmds
# ------------------------------------------------------------------------------
# 2021-2024 data
oyster_biometry_2021_2024 <- read_excel(
old_file,
sheet = "oyster_biometry"
) %>%
clean_names()
tax_counts_2021_2024 <- read_excel(
old_file,
sheet = "taxonomic_group_counts"
) %>%
clean_names()
# 2025 data
oyster_biometry_2025 <- read_excel(
new_file,
sheet = "oyster_biometry"
) %>%
clean_names()
samples_2025 <- read_excel(
new_file,
sheet = "samples"
) %>%
clean_names()
# Taxa recorded in the 2021-2024 taxonomic count sheet.
# These exclude live/dead oysters because the goal here is associated biodiversity.
taxa_cols <- c(
"algae", "barnacle", "mangrove_crab", "mud_crab", "swamp_eel",
"sponge", "uca_crab", "blue_mussel", "polychaete", "mussel",
"sea_snail", "swimming_crab", "crab", "hermit_crab",
"mangrove_tree_crab", "snapping_shrimp", "fish"
)
taxa_cols <- taxa_cols[taxa_cols %in% names(tax_counts_2021_2024)]
# ------------------------------------------------------------------------------
# 2. Prepare density and live/dead oyster data
# ------------------------------------------------------------------------------
clean_density_fields <- function(df, source_label) {
# Some sheets do not have these columns; adding them avoids failures.
for (col in c("sampling_year", "sampling_month")) {
if (!col %in% names(df)) df[[col]] <- NA
}
df %>%
mutate(
data_source = source_label,
oyster_reef_name = normalize_reef(oyster_reef_name),
seasonal_period = normalize_season(seasonal_period),
sampling_date = to_sampling_date(sampling_day),
sampling_year_from_col = as.integer(parse_numeric_safely(sampling_year)),
sampling_year = coalesce(sampling_year_from_col, year(sampling_date)),
sampling_month = coalesce(
as.character(sampling_month),
as.character(month(sampling_date, label = TRUE, abbr = FALSE))
),
sampling_month = factor(sampling_month, levels = month.name),
stratum_no = parse_numeric_safely(stratum_no),
quadrat_no = parse_numeric_safely(quadrat_no),
live_oysters = parse_numeric_safely(live_oysters),
dead_oysters = parse_numeric_safely(dead_oysters),
live_oysters = replace_na(live_oysters, 0),
dead_oysters = replace_na(dead_oysters, 0),
total_oysters = live_oysters + dead_oysters,
live_prop_quadrat = if_else(total_oysters > 0, live_oysters / total_oysters, NA_real_),
dead_prop_quadrat = if_else(total_oysters > 0, dead_oysters / total_oysters, NA_real_),
total_density = total_oysters / quadrat_area_m2,
live_density = live_oysters / quadrat_area_m2,
dead_density = dead_oysters / quadrat_area_m2
) %>%
filter(!is.na(oyster_reef_name), !is.na(sampling_year)) %>%
select(
data_source,
country, ma_name, community, waterbody_name, oyster_reef_name,
sampling_date, sampling_year, sampling_month, seasonal_period, sampling_time,
stratum_no, quadrat_no,
live_oysters, dead_oysters, total_oysters,
live_prop_quadrat, dead_prop_quadrat,
total_density, live_density, dead_density
)
}
density_2021_2024 <- clean_density_fields(
tax_counts_2021_2024,
source_label = "2021-2024 taxonomic_group_counts"
)
density_2025 <- clean_density_fields(
samples_2025,
source_label = "2025 samples"
)
density_data <- bind_rows(density_2021_2024, density_2025) %>%
mutate(
sampling_year = as.integer(sampling_year),
oyster_reef_name = fct_reorder(oyster_reef_name, sampling_year, .fun = min, .na_rm = TRUE)
)
# Helper to summarise means and 95% confidence intervals
summarise_density <- function(data, group_vars) {
data %>%
group_by(across(all_of(group_vars))) %>%
summarise(
n_quadrats = n(),
n_strata = n_distinct(stratum_no, na.rm = TRUE),
mean_density = mean(total_density, na.rm = TRUE),
sd_density = sd(total_density, na.rm = TRUE),
se_density = if_else(n_quadrats > 1, sd_density / sqrt(n_quadrats), NA_real_),
t_crit = if_else(n_quadrats > 1, qt(0.975, df = n_quadrats - 1), NA_real_),
ci_lower_density = pmax(0, mean_density - t_crit * se_density),
ci_upper_density = mean_density + t_crit * se_density,
mean_live_density = mean(live_density, na.rm = TRUE),
mean_dead_density = mean(dead_density, na.rm = TRUE),
total_live_oysters = sum(live_oysters, na.rm = TRUE),
total_dead_oysters = sum(dead_oysters, na.rm = TRUE),
total_oysters_counted = total_live_oysters + total_dead_oysters,
.groups = "drop"
) %>%
mutate(
interpretation = case_when(
n_quadrats >= min_quadrats_seasonal ~ "Good for reef-level/seasonal interpretation",
n_quadrats >= min_quadrats_main ~ "Usable descriptively; interpret seasonality cautiously",
TRUE ~ "Low sample size; avoid strong interpretation"
)
)
}
density_reef_year <- summarise_density(
density_data,
c("oyster_reef_name", "sampling_year")
)
density_reef_year_season <- summarise_density(
density_data,
c("oyster_reef_name", "sampling_year", "seasonal_period")
) %>%
filter(!is.na(seasonal_period))
# Live/dead proportions calculated as pooled oyster-level proportions within each group.
# This is usually more intuitive for presentation than averaging quadrat-level proportions.
summarise_live_dead <- function(data, group_vars) {
data %>%
group_by(across(all_of(group_vars))) %>%
summarise(
n_quadrats = n(),
live_oysters = sum(live_oysters, na.rm = TRUE),
dead_oysters = sum(dead_oysters, na.rm = TRUE),
total_oysters = live_oysters + dead_oysters,
.groups = "drop"
) %>%
mutate(
live_prop = if_else(total_oysters > 0, live_oysters / total_oysters, NA_real_),
dead_prop = if_else(total_oysters > 0, dead_oysters / total_oysters, NA_real_),
live_ci_low = wilson_ci_low(live_oysters, total_oysters),
live_ci_high = wilson_ci_high(live_oysters, total_oysters),
interpretation = case_when(
n_quadrats >= min_quadrats_seasonal ~ "Good for reef-level/seasonal interpretation",
n_quadrats >= min_quadrats_main ~ "Usable descriptively; interpret seasonality cautiously",
TRUE ~ "Low sample size; avoid strong interpretation"
)
)
}
live_dead_reef_year <- summarise_live_dead(
density_data,
c("oyster_reef_name", "sampling_year")
)
live_dead_reef_year_season <- summarise_live_dead(
density_data,
c("oyster_reef_name", "sampling_year", "seasonal_period")
) %>%
filter(!is.na(seasonal_period))
# Seasonal statistical comparisons.
# These are optional and should be used as supporting evidence, not as the main message.
seasonal_density_tests <- density_data %>%
filter(seasonal_period %in% c("Dry", "Rainy")) %>%
group_by(oyster_reef_name, sampling_year) %>%
group_modify(~ {
n_dry <- sum(.x$seasonal_period == "Dry", na.rm = TRUE)
n_rainy <- sum(.x$seasonal_period == "Rainy", na.rm = TRUE)
ok <- n_dry >= min_quadrats_seasonal && n_rainy >= min_quadrats_seasonal
tibble(
n_dry = n_dry,
n_rainy = n_rainy,
test_used = if_else(ok, "Wilcoxon rank-sum", "Not tested: low/incomplete seasonal sample"),
p_value = if (ok) {
wilcox.test(total_density ~ seasonal_period, data = .x)$p.value
} else {
NA_real_
}
)
}) %>%
ungroup()
seasonal_density_contrast <- density_reef_year_season %>%
filter(seasonal_period %in% c("Dry", "Rainy")) %>%
select(oyster_reef_name, sampling_year, seasonal_period, n_quadrats, mean_density) %>%
pivot_wider(
names_from = seasonal_period,
values_from = c(mean_density, n_quadrats),
names_sep = "_"
) %>%
mutate(
density_difference_rainy_minus_dry = mean_density_Rainy - mean_density_Dry,
density_pct_difference_rainy_vs_dry = if_else(
!is.na(mean_density_Dry) & mean_density_Dry > 0,
100 * density_difference_rainy_minus_dry / mean_density_Dry,
NA_real_
),
density_pct_difference_label = if_else(
is.na(density_pct_difference_rainy_vs_dry),
"",
paste0(round(density_pct_difference_rainy_vs_dry), "%")
),
seasonal_contrast_interpretation = case_when(
n_quadrats_Dry >= min_quadrats_seasonal &
n_quadrats_Rainy >= min_quadrats_seasonal ~
"Seasonal contrast has adequate quadrat coverage",
TRUE ~ "Seasonal contrast should be interpreted cautiously"
)
) %>%
left_join(seasonal_density_tests, by = c("oyster_reef_name", "sampling_year"))
# ------------------------------------------------------------------------------
# 3. Prepare oyster size data for condition relationships
# ------------------------------------------------------------------------------
clean_biometry_fields <- function(df, source_label) {
for (col in c("sampling_year", "sampling_month", "commercial_size_class_mm")) {
if (!col %in% names(df)) df[[col]] <- NA
}
df %>%
mutate(
data_source = source_label,
oyster_reef_name = normalize_reef(oyster_reef_name),
seasonal_period = normalize_season(seasonal_period),
sampling_date = to_sampling_date(sampling_day),
sampling_year_from_col = as.integer(parse_numeric_safely(sampling_year)),
sampling_year = coalesce(sampling_year_from_col, year(sampling_date)),
height_mm = parse_numeric_safely(height_mm),
length_mm = parse_numeric_safely(length_mm),
size_class_length = case_when(
!is.na(length_mm) & length_mm <= 29 ~ "Seed (≤29 mm)",
!is.na(length_mm) & length_mm >= 30 & length_mm <= 59 ~ "Juvenile (30–59 mm)",
!is.na(length_mm) & length_mm >= 60 ~ "Commercial-size (≥60 mm)",
TRUE ~ NA_character_
),
size_class_length = factor(
size_class_length,
levels = c("Seed (≤29 mm)", "Juvenile (30–59 mm)", "Commercial-size (≥60 mm)")
)
) %>%
filter(!is.na(oyster_reef_name), !is.na(sampling_year)) %>%
select(
data_source, country, ma_name, community, waterbody_name, oyster_reef_name,
sampling_date, sampling_year, sampling_month, seasonal_period, sampling_time,
stratum_no, quadrat_no, height_mm, length_mm,
commercial_size_class_mm, size_class_length
)
}
all_oyster_biometry <- bind_rows(
clean_biometry_fields(oyster_biometry_2021_2024, "2021-2024 oyster_biometry"),
clean_biometry_fields(oyster_biometry_2025, "2025 oyster_biometry")
)
size_summary_reef_year <- all_oyster_biometry %>%
group_by(oyster_reef_name, sampling_year) %>%
summarise(
n_measured_oysters = n(),
mean_height_mm = mean(height_mm, na.rm = TRUE),
sd_height_mm = sd(height_mm, na.rm = TRUE),
mean_length_mm = mean(length_mm, na.rm = TRUE),
sd_length_mm = sd(length_mm, na.rm = TRUE),
prop_commercial_size = mean(length_mm >= 60, na.rm = TRUE),
.groups = "drop"
)
size_structure_reef_year <- all_oyster_biometry %>%
filter(!is.na(size_class_length)) %>%
count(oyster_reef_name, sampling_year, size_class_length, name = "n_oysters") %>%
group_by(oyster_reef_name, sampling_year) %>%
mutate(prop_oysters = n_oysters / sum(n_oysters)) %>%
ungroup()
# ------------------------------------------------------------------------------
# 4. Prepare associated biodiversity / taxonomic richness data
# ------------------------------------------------------------------------------
# Note: associated taxa richness is currently only calculated from the 2021-2024
# taxonomic_group_counts sheet because the 2025 workbook does not include the same
# associated taxa columns in the samples sheet.
tax_counts_clean <- tax_counts_2021_2024 %>%
mutate(
oyster_reef_name = normalize_reef(oyster_reef_name),
seasonal_period = normalize_season(seasonal_period),
sampling_date = to_sampling_date(sampling_day),
sampling_year = year(sampling_date),
sampling_month = factor(month(sampling_date, label = TRUE, abbr = FALSE), levels = month.name),
stratum_no = parse_numeric_safely(stratum_no),
quadrat_no = parse_numeric_safely(quadrat_no),
across(any_of(c("live_oysters", "dead_oysters", taxa_cols)), parse_numeric_safely)
) %>%
filter(!is.na(oyster_reef_name), !is.na(sampling_year)) %>%
mutate(
across(any_of(taxa_cols), ~ replace_na(.x, 0)),
associated_taxa_richness_quadrat = rowSums(across(all_of(taxa_cols), ~ .x > 0), na.rm = TRUE),
associated_taxa_abundance_quadrat = rowSums(across(all_of(taxa_cols)), na.rm = TRUE)
)
taxa_long <- tax_counts_clean %>%
select(
oyster_reef_name, sampling_year, sampling_month, seasonal_period,
stratum_no, quadrat_no, all_of(taxa_cols)
) %>%
pivot_longer(
cols = all_of(taxa_cols),
names_to = "taxon",
values_to = "taxon_count"
) %>%
mutate(
taxon_count = replace_na(taxon_count, 0),
taxon_present = taxon_count > 0,
taxon_label = taxon %>% str_replace_all("_", " ") %>% str_to_sentence()
)
biodiversity_summary_reef_year_base <- tax_counts_clean %>%
group_by(oyster_reef_name, sampling_year) %>%
summarise(
n_quadrats_biodiversity = n(),
mean_associated_taxa_richness_per_quadrat = mean(associated_taxa_richness_quadrat, na.rm = TRUE),
sd_associated_taxa_richness_per_quadrat = sd(associated_taxa_richness_quadrat, na.rm = TRUE),
mean_associated_taxa_abundance_per_quadrat = mean(associated_taxa_abundance_quadrat, na.rm = TRUE),
total_associated_taxa_abundance = sum(associated_taxa_abundance_quadrat, na.rm = TRUE),
observed_associated_taxa_richness = sum(map_lgl(across(all_of(taxa_cols)), ~ any(replace_na(.x, 0) > 0))),
.groups = "drop"
)
taxa_totals_reef_year <- taxa_long %>%
group_by(oyster_reef_name, sampling_year, taxon, taxon_label) %>%
summarise(
total_count = sum(taxon_count, na.rm = TRUE),
occurrence_quadrats = sum(taxon_present, na.rm = TRUE),
n_quadrats = n(),
occurrence_pct = occurrence_quadrats / n_quadrats,
.groups = "drop"
)
shannon_reef_year <- taxa_totals_reef_year %>%
group_by(oyster_reef_name, sampling_year) %>%
summarise(
shannon_associated_taxa = shannon_index(total_count),
.groups = "drop"
)
biodiversity_summary_reef_year <- biodiversity_summary_reef_year_base %>%
left_join(shannon_reef_year, by = c("oyster_reef_name", "sampling_year"))
# Combined reef-year dataset for exploratory relationship plots.
# These relationships should be described as exploratory correlations, not causality.
condition_biodiversity <- biodiversity_summary_reef_year %>%
left_join(
density_reef_year %>%
select(
oyster_reef_name, sampling_year, n_quadrats_density = n_quadrats,
mean_density, ci_lower_density, ci_upper_density,
mean_live_density, mean_dead_density, total_oysters_counted
),
by = c("oyster_reef_name", "sampling_year")
) %>%
left_join(
live_dead_reef_year %>%
select(oyster_reef_name, sampling_year, live_prop, dead_prop),
by = c("oyster_reef_name", "sampling_year")
) %>%
left_join(
size_summary_reef_year,
by = c("oyster_reef_name", "sampling_year")
)
relationship_correlations <- condition_biodiversity %>%
summarise(
n_reef_years = n(),
spearman_observed_richness_vs_density = safe_cor(observed_associated_taxa_richness, mean_density),
spearman_mean_quadrat_richness_vs_density = safe_cor(mean_associated_taxa_richness_per_quadrat, mean_density),
spearman_observed_richness_vs_live_prop = safe_cor(observed_associated_taxa_richness, live_prop),
spearman_observed_richness_vs_mean_length = safe_cor(observed_associated_taxa_richness, mean_length_mm),
spearman_shannon_vs_density = safe_cor(shannon_associated_taxa, mean_density),
spearman_shannon_vs_live_prop = safe_cor(shannon_associated_taxa, live_prop)
)
# ------------------------------------------------------------------------------
# 5. Tables for review and/or presentation backup
# ------------------------------------------------------------------------------
sampling_interpretability_table <- density_reef_year_season %>%
select(
oyster_reef_name, sampling_year, seasonal_period, n_quadrats,
mean_density, ci_lower_density, ci_upper_density, interpretation
) %>%
arrange(oyster_reef_name, sampling_year, seasonal_period)
write_csv(density_reef_year, file.path(table_dir, "density_by_reef_year.csv"))
write_csv(density_reef_year_season, file.path(table_dir, "density_by_reef_year_season.csv"))
write_csv(live_dead_reef_year, file.path(table_dir, "live_dead_proportions_by_reef_year.csv"))
write_csv(live_dead_reef_year_season, file.path(table_dir, "live_dead_proportions_by_reef_year_season.csv"))
write_csv(seasonal_density_contrast, file.path(table_dir, "seasonal_density_contrast_rainy_vs_dry.csv"))
write_csv(sampling_interpretability_table, file.path(table_dir, "seasonal_interpretability_table.csv"))
write_csv(size_summary_reef_year, file.path(table_dir, "size_summary_by_reef_year.csv"))
write_csv(size_structure_reef_year, file.path(table_dir, "size_structure_by_reef_year.csv"))
write_csv(biodiversity_summary_reef_year, file.path(table_dir, "biodiversity_summary_by_reef_year_2021_2024.csv"))
write_csv(taxa_totals_reef_year, file.path(table_dir, "taxa_occurrence_and_counts_by_reef_year_2021_2024.csv"))
write_csv(condition_biodiversity, file.path(table_dir, "condition_biodiversity_relationship_dataset_2021_2024.csv"))
write_csv(relationship_correlations, file.path(table_dir, "condition_biodiversity_spearman_correlations.csv"))
# ------------------------------------------------------------------------------
# 6. Plot set A: Sampling effort and density
# ------------------------------------------------------------------------------
p_sampling_effort <- density_data %>%
count(oyster_reef_name, sampling_year, seasonal_period, name = "n_quadrats") %>%
filter(!is.na(seasonal_period)) %>%
ggplot(aes(x = factor(sampling_year), y = fct_rev(oyster_reef_name), fill = n_quadrats)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = n_quadrats), size = 3) +
facet_wrap(~ seasonal_period) +
scale_fill_viridis_c(option = "C", direction = -1) +
labs(
title = "Sampling effort by reef, year, and season",
subtitle = "Number of quadrats available for density and live/dead analyses",
x = NULL,
y = NULL,
fill = "Quadrats",
caption = paste0(
"Rule of thumb used here: ≥", min_quadrats_seasonal,
" quadrats = stronger seasonal interpretation; ", min_quadrats_main,
"–", min_quadrats_seasonal - 1, " = descriptive/cautious."
)
) +
theme_oyster()
p_density_year <- density_reef_year %>%
ggplot(aes(x = sampling_year, y = mean_density, group = oyster_reef_name)) +
geom_ribbon(
aes(ymin = ci_lower_density, ymax = ci_upper_density),
alpha = 0.12
) +
geom_line(linewidth = 0.6) +
geom_point(aes(size = n_quadrats), alpha = 0.85) +
facet_wrap(~ oyster_reef_name, scales = "free_y") +
scale_x_continuous(breaks = sort(unique(density_reef_year$sampling_year))) +
scale_size_continuous(range = c(1.5, 5)) +
labs(
title = "Mean oyster density by reef and year",
subtitle = paste0("Mean density with 95% confidence intervals; unit = ", density_unit_label),
x = NULL,
y = density_unit_label,
size = "Quadrats",
caption = "Density is calculated as live + dead oysters. If the quadrat area is updated, this plot automatically switches to density per m²."
) +
theme_oyster()
p_density_heatmap <- density_reef_year %>%
ggplot(aes(x = factor(sampling_year), y = fct_rev(oyster_reef_name), fill = mean_density)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = round(mean_density, 1)), size = 3) +
scale_fill_viridis_c(option = "mako", direction = -1) +
labs(
title = "Oyster density heatmap",
subtitle = "Reef-year comparison",
x = NULL,
y = NULL,
fill = density_unit_label
) +
theme_oyster()
save_plot(p_sampling_effort, "01_sampling_effort_reef_year_season.png", width = 12, height = 7)
save_plot(p_density_year, "02_mean_density_by_reef_year_ci.png", width = 14, height = 9)
save_plot(p_density_heatmap, "03_density_heatmap_by_reef_year.png", width = 10, height = 7)
# ------------------------------------------------------------------------------
# 7. Plot set B: Seasonal comparisons
# ------------------------------------------------------------------------------
p_density_season <- density_reef_year_season %>%
ggplot(aes(x = factor(sampling_year), y = mean_density, color = seasonal_period)) +
geom_errorbar(
aes(ymin = ci_lower_density, ymax = ci_upper_density),
width = 0.15,
position = position_dodge(width = 0.5),
alpha = 0.75
) +
geom_point(
aes(shape = interpretation, size = n_quadrats),
position = position_dodge(width = 0.5),
alpha = 0.9
) +
facet_wrap(~ oyster_reef_name, scales = "free_y") +
labs(
title = "Seasonal oyster density by reef and year",
subtitle = "Rainy vs dry comparisons are flagged by sampling coverage",
x = NULL,
y = density_unit_label,
color = "Season",
shape = "Interpretation flag",
size = "Quadrats"
) +
theme_oyster(base_size = 11)
p_seasonal_contrast <- seasonal_density_contrast %>%
ggplot(aes(x = factor(sampling_year), y = fct_rev(oyster_reef_name), fill = density_pct_difference_rainy_vs_dry)) +
geom_tile(color = "white", linewidth = 0.4) +
geom_text(aes(label = if_else(is.na(density_pct_difference_rainy_vs_dry), "", paste0(round(density_pct_difference_rainy_vs_dry), "%"))), size = 3) +
scale_fill_gradient2(labels = label_percent(scale = 1), na.value = "grey90") +
labs(
title = "Seasonal density contrast: rainy vs dry",
subtitle = "Positive values mean higher density in the rainy period; negative values mean higher density in the dry period",
x = NULL,
y = NULL,
fill = "Rainy vs dry\n% difference",
caption = "Blank cells indicate that one of the seasons was not available for that reef-year."
) +
theme_oyster()
save_plot(p_density_season, "04_seasonal_density_by_reef_year_ci.png", width = 14, height = 9)
save_plot(p_seasonal_contrast, "05_seasonal_density_contrast_heatmap.png", width = 11, height = 7)
# ------------------------------------------------------------------------------
# 8. Plot set C: Live/dead oyster proportions
# ------------------------------------------------------------------------------
live_dead_counts_long <- live_dead_reef_year %>%
select(oyster_reef_name, sampling_year, live_oysters, dead_oysters, n_quadrats) %>%
pivot_longer(
cols = c(live_oysters, dead_oysters),
names_to = "status",
values_to = "oyster_count"
) %>%
mutate(
status = recode(status, live_oysters = "Live oysters", dead_oysters = "Dead oysters"),
status = factor(status, levels = c("Live oysters", "Dead oysters"))
)
p_live_dead_stack <- live_dead_counts_long %>%
ggplot(aes(x = factor(sampling_year), y = oyster_count, fill = status)) +
geom_col(position = "fill", width = 0.8) +
facet_wrap(~ oyster_reef_name) +
scale_y_continuous(labels = label_percent()) +
scale_fill_manual(values = c("Live oysters" = "#2a9d8f", "Dead oysters" = "#e76f51")) +
labs(
title = "Proportion of live and dead oysters by reef and year",
subtitle = "Pooled oyster-level proportion from quadrat counts",
x = NULL,
y = "Proportion of counted oysters",
fill = NULL
) +
theme_oyster()
p_live_prop_ci <- live_dead_reef_year %>%
ggplot(aes(x = sampling_year, y = live_prop, group = oyster_reef_name)) +
geom_ribbon(aes(ymin = live_ci_low, ymax = live_ci_high), alpha = 0.12) +
geom_line(linewidth = 0.6) +
geom_point(aes(size = n_quadrats), alpha = 0.9) +
facet_wrap(~ oyster_reef_name) +
scale_x_continuous(breaks = sort(unique(live_dead_reef_year$sampling_year))) +
scale_y_continuous(labels = label_percent(), limits = c(0, 1)) +
labs(
title = "Live oyster proportion by reef and year",
subtitle = "Points show live proportion; ribbons show 95% confidence intervals",
x = NULL,
y = "Live oyster proportion",
size = "Quadrats"
) +
theme_oyster()
p_live_dead_season <- live_dead_reef_year_season %>%
select(oyster_reef_name, sampling_year, seasonal_period, live_oysters, dead_oysters, n_quadrats) %>%
pivot_longer(
cols = c(live_oysters, dead_oysters),
names_to = "status",
values_to = "oyster_count"
) %>%
mutate(
status = recode(status, live_oysters = "Live oysters", dead_oysters = "Dead oysters"),
status = factor(status, levels = c("Live oysters", "Dead oysters"))
) %>%
ggplot(aes(x = seasonal_period, y = oyster_count, fill = status)) +
geom_col(position = "fill") +
facet_grid(oyster_reef_name ~ sampling_year) +
scale_y_continuous(labels = label_percent()) +
scale_fill_manual(values = c("Live oysters" = "#2a9d8f", "Dead oysters" = "#e76f51")) +
labs(
title = "Seasonal live/dead oyster composition",
subtitle = "Use cautiously where seasonal quadrat coverage is low",
x = NULL,
y = "Proportion of counted oysters",
fill = NULL
) +
theme_oyster(base_size = 10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
save_plot(p_live_dead_stack, "06_live_dead_proportions_by_reef_year.png", width = 13, height = 8)
save_plot(p_live_prop_ci, "07_live_oyster_proportion_ci_by_reef_year.png", width = 13, height = 8)
save_plot(p_live_dead_season, "08_seasonal_live_dead_proportions.png", width = 15, height = 12)
# ------------------------------------------------------------------------------
# 9. Plot set D: Associated biodiversity and taxonomic richness
# ------------------------------------------------------------------------------
p_richness_year <- biodiversity_summary_reef_year %>%
ggplot(aes(x = sampling_year, y = observed_associated_taxa_richness, group = oyster_reef_name)) +
geom_line(linewidth = 0.6) +
geom_point(aes(size = n_quadrats_biodiversity), alpha = 0.9) +
facet_wrap(~ oyster_reef_name) +
scale_x_continuous(breaks = sort(unique(biodiversity_summary_reef_year$sampling_year))) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(
title = "Observed associated taxa richness by reef and year",
subtitle = "Number of associated taxa observed at least once in the reef-year sample",
x = NULL,
y = "Observed associated taxa richness",
size = "Quadrats"
) +
theme_oyster()
p_mean_richness_quadrat <- biodiversity_summary_reef_year %>%
ggplot(aes(x = sampling_year, y = mean_associated_taxa_richness_per_quadrat, group = oyster_reef_name)) +
geom_line(linewidth = 0.6) +
geom_point(aes(size = n_quadrats_biodiversity), alpha = 0.9) +
facet_wrap(~ oyster_reef_name) +
scale_x_continuous(breaks = sort(unique(biodiversity_summary_reef_year$sampling_year))) +
labs(
title = "Mean associated taxa richness per quadrat",
subtitle = "Average number of associated taxa present per quadrat",
x = NULL,
y = "Mean taxa richness per quadrat",
size = "Quadrats"
) +
theme_oyster()
p_taxa_occurrence_heatmap <- taxa_totals_reef_year %>%
group_by(oyster_reef_name, taxon_label) %>%
summarise(
mean_occurrence_pct = mean(occurrence_pct, na.rm = TRUE),
total_count = sum(total_count, na.rm = TRUE),
.groups = "drop"
) %>%
filter(total_count > 0) %>%
ggplot(aes(x = taxon_label, y = fct_rev(oyster_reef_name), fill = mean_occurrence_pct)) +
geom_tile(color = "white", linewidth = 0.35) +
scale_fill_viridis_c(labels = label_percent(), option = "C") +
labs(
title = "Associated taxa occurrence by reef",
subtitle = "Average proportion of quadrats where each associated taxon was observed across 2021-2024",
x = NULL,
y = NULL,
fill = "Occurrence"
) +
theme_oyster(base_size = 11) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_taxa_composition <- taxa_totals_reef_year %>%
filter(total_count > 0) %>%
ggplot(aes(x = factor(sampling_year), y = total_count, fill = taxon_label)) +
geom_col(position = "fill") +
facet_wrap(~ oyster_reef_name) +
scale_y_continuous(labels = label_percent()) +
labs(
title = "Associated biodiversity composition by reef and year",
subtitle = "Relative composition of associated taxa counts",
x = NULL,
y = "Share of associated taxa counts",
fill = "Taxon"
) +
theme_oyster(base_size = 10)
save_plot(p_richness_year, "09_observed_associated_taxa_richness_by_reef_year.png", width = 13, height = 8)
save_plot(p_mean_richness_quadrat, "10_mean_taxa_richness_per_quadrat_by_reef_year.png", width = 13, height = 8)
save_plot(p_taxa_occurrence_heatmap, "11_taxa_occurrence_heatmap_by_reef.png", width = 13, height = 7)
save_plot(p_taxa_composition, "12_taxa_composition_by_reef_year.png", width = 14, height = 9)
# ------------------------------------------------------------------------------
# 10. Plot set E: Exploratory biodiversity-condition relationships
# ------------------------------------------------------------------------------
p_relationship_density <- condition_biodiversity %>%
ggplot(aes(x = mean_density, y = observed_associated_taxa_richness)) +
geom_point(aes(size = n_quadrats_biodiversity, color = factor(sampling_year)), alpha = 0.85) +
geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
geom_text(aes(label = oyster_reef_name), check_overlap = TRUE, vjust = -0.7, size = 3) +
labs(
title = "Associated taxa richness vs oyster density",
subtitle = "Exploratory reef-year relationship; use to guide discussion, not causal inference",
x = density_unit_label,
y = "Observed associated taxa richness",
color = "Year",
size = "Biodiversity quadrats"
) +
theme_oyster()
p_relationship_live_prop <- condition_biodiversity %>%
ggplot(aes(x = live_prop, y = observed_associated_taxa_richness)) +
geom_point(aes(size = n_quadrats_biodiversity, color = factor(sampling_year)), alpha = 0.85) +
geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
geom_text(aes(label = oyster_reef_name), check_overlap = TRUE, vjust = -0.7, size = 3) +
scale_x_continuous(labels = label_percent()) +
labs(
title = "Associated taxa richness vs live oyster proportion",
subtitle = "Exploratory reef-year relationship between associated biodiversity and oyster condition",
x = "Live oyster proportion",
y = "Observed associated taxa richness",
color = "Year",
size = "Biodiversity quadrats"
) +
theme_oyster()
p_relationship_size <- condition_biodiversity %>%
ggplot(aes(x = mean_length_mm, y = observed_associated_taxa_richness)) +
geom_point(aes(size = n_quadrats_biodiversity, color = factor(sampling_year)), alpha = 0.85) +
geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
geom_text(aes(label = oyster_reef_name), check_overlap = TRUE, vjust = -0.7, size = 3) +
labs(
title = "Associated taxa richness vs mean oyster length",
subtitle = "Exploratory reef-year relationship between associated biodiversity and oyster size structure",
x = "Mean oyster length (mm)",
y = "Observed associated taxa richness",
color = "Year",
size = "Biodiversity quadrats"
) +
theme_oyster()
p_relationship_shannon_density <- condition_biodiversity %>%
ggplot(aes(x = mean_density, y = shannon_associated_taxa)) +
geom_point(aes(size = n_quadrats_biodiversity, color = factor(sampling_year)), alpha = 0.85) +
geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
geom_text(aes(label = oyster_reef_name), check_overlap = TRUE, vjust = -0.7, size = 3) +
labs(
title = "Associated taxa diversity vs oyster density",
subtitle = "Shannon index calculated from associated taxa counts by reef-year",
x = density_unit_label,
y = "Shannon associated taxa index",
color = "Year",
size = "Biodiversity quadrats"
) +
theme_oyster()
save_plot(p_relationship_density, "13_relationship_richness_vs_oyster_density.png", width = 10, height = 7)
save_plot(p_relationship_live_prop, "14_relationship_richness_vs_live_oyster_proportion.png", width = 10, height = 7)
save_plot(p_relationship_size, "15_relationship_richness_vs_mean_oyster_length.png", width = 10, height = 7)
save_plot(p_relationship_shannon_density, "16_relationship_shannon_vs_oyster_density.png", width = 10, height = 7)
# ------------------------------------------------------------------------------
# 11. Plot set F: Oyster size structure as reef-condition context
# ------------------------------------------------------------------------------
p_size_structure <- size_structure_reef_year %>%
ggplot(aes(x = factor(sampling_year), y = prop_oysters, fill = size_class_length)) +
geom_col(width = 0.8) +
facet_wrap(~ oyster_reef_name) +
scale_y_continuous(labels = label_percent()) +
scale_fill_manual(
values = c(
"Seed (≤29 mm)" = "#8ecae6",
"Juvenile (30–59 mm)" = "#219ebc",
"Commercial-size (≥60 mm)" = "#023047"
)
) +
labs(
title = "Oyster size structure by reef and year",
subtitle = "Length-based size classes using thresholds from the 2025 commercial-size classification",
x = NULL,
y = "Proportion of measured oysters",
fill = "Size class"
) +
theme_oyster()
p_mean_length_density <- size_summary_reef_year %>%
left_join(
density_reef_year %>% select(oyster_reef_name, sampling_year, mean_density, n_quadrats),
by = c("oyster_reef_name", "sampling_year")
) %>%
ggplot(aes(x = mean_density, y = mean_length_mm)) +
geom_point(aes(size = n_quadrats, color = factor(sampling_year)), alpha = 0.85) +
geom_smooth(method = "lm", se = FALSE, color = "grey30", linewidth = 0.7) +
geom_text(aes(label = oyster_reef_name), check_overlap = TRUE, vjust = -0.7, size = 3) +
labs(
title = "Mean oyster length vs density",
subtitle = "Useful for discussing whether denser reefs also show different size structure",
x = density_unit_label,
y = "Mean oyster length (mm)",
color = "Year",
size = "Quadrats"
) +
theme_oyster()
save_plot(p_size_structure, "17_oyster_size_structure_by_reef_year.png", width = 13, height = 8)
save_plot(p_mean_length_density, "18_relationship_mean_length_vs_density.png", width = 10, height = 7)
# ------------------------------------------------------------------------------
# 12. Optional tables formatted for Rmd display
# ------------------------------------------------------------------------------
# If this code is pasted into an Rmd, these example tables are useful for quick review.
# Uncomment the chunks below inside an Rmd if you want pretty HTML tables.
# density_reef_year %>%
# select(
# Reef = oyster_reef_name,
# Year = sampling_year,
# `Quadrats` = n_quadrats,
# `Mean density` = mean_density,
# `Lower 95% CI` = ci_lower_density,
# `Upper 95% CI` = ci_upper_density,
# `Interpretation` = interpretation
# ) %>%
# mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
# kbl(caption = paste0("Mean oyster density by reef and year (", density_unit_label, ")."), booktabs = TRUE) %>%
# kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
# live_dead_reef_year %>%
# select(
# Reef = oyster_reef_name,
# Year = sampling_year,
# `Quadrats` = n_quadrats,
# `Live oysters` = live_oysters,
# `Dead oysters` = dead_oysters,
# `Live proportion` = live_prop,
# `Lower 95% CI` = live_ci_low,
# `Upper 95% CI` = live_ci_high
# ) %>%
# mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
# kbl(caption = "Live oyster proportion by reef and year.", booktabs = TRUE) %>%
# kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
# ------------------------------------------------------------------------------
# 13. Plot index / what to use in the presentation
# ------------------------------------------------------------------------------
plot_index <- tribble(
~file, ~best_use,
"01_sampling_effort_reef_year_season.png", "Start here: shows where seasonal interpretation is strong vs weak.",
"02_mean_density_by_reef_year_ci.png", "Main density result by reef and year with uncertainty.",
"03_density_heatmap_by_reef_year.png", "Compact slide-friendly density comparison.",
"04_seasonal_density_by_reef_year_ci.png", "Seasonal density comparison, best for technical discussion.",
"05_seasonal_density_contrast_heatmap.png", "Simple rainy-vs-dry contrast; use only with sample-size caveat.",
"06_live_dead_proportions_by_reef_year.png", "Main live/dead composition plot by reef and year.",
"07_live_oyster_proportion_ci_by_reef_year.png", "Live oyster condition trend with uncertainty.",
"08_seasonal_live_dead_proportions.png", "Seasonal live/dead composition; use cautiously.",
"09_observed_associated_taxa_richness_by_reef_year.png", "Observed associated biodiversity by reef-year, 2021-2024 only.",
"10_mean_taxa_richness_per_quadrat_by_reef_year.png", "Average richness per quadrat, less sensitive to total sampling effort.",
"11_taxa_occurrence_heatmap_by_reef.png", "Best biodiversity slide for community discussion: which taxa appear where.",
"12_taxa_composition_by_reef_year.png", "Composition of associated taxa counts by reef-year.",
"13_relationship_richness_vs_oyster_density.png", "Exploratory relationship between biodiversity and reef density.",
"14_relationship_richness_vs_live_oyster_proportion.png", "Exploratory relationship between biodiversity and live oyster proportion.",
"15_relationship_richness_vs_mean_oyster_length.png", "Exploratory relationship between biodiversity and size structure.",
"16_relationship_shannon_vs_oyster_density.png", "Alternative biodiversity metric: Shannon diversity vs density.",
"17_oyster_size_structure_by_reef_year.png", "Size-structure context for reef condition.",
"18_relationship_mean_length_vs_density.png", "Density-size relationship for discussion of reef condition."
)
write_csv(plot_index, file.path(output_dir, "plot_index.csv"))
message("Done. Outputs saved to: ", output_dir)
message("Density unit used: ", density_unit_label)
message("Remember to confirm quadrat_area_m2 before presenting density per m².")
1. Data coverage and interpretation flags
Before interpreting seasonal patterns, maybe it is useful to know where sample coverage is strong enough to support discussion and where results should remain descriptive.
| Reef | Year | Season | Number of quadrats | Mean oyster density | Lower 95% CI | Upper 95% CI | Interpretation flag |
|---|---|---|---|---|---|---|---|
| Água Boa | 2021 | Dry | 16 | 33.69 | 21.80 | 45.58 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2022 | Dry | 15 | 12.93 | 6.57 | 19.30 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2022 | Rainy | 15 | 21.00 | 10.14 | 31.86 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2023 | Dry | 15 | 14.00 | 7.67 | 20.33 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2023 | Rainy | 15 | 9.60 | 4.09 | 15.11 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2024 | Dry | 15 | 11.20 | 5.51 | 16.89 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2024 | Rainy | 15 | 16.27 | 7.95 | 24.59 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2025 | Dry | 20 | 5.20 | 0.82 | 9.58 | Good for reef-level/seasonal interpretation |
| Água Boa | 2025 | Rainy | 30 | 8.33 | 3.08 | 13.58 | Good for reef-level/seasonal interpretation |
| Aquavila | 2021 | Dry | 40 | 45.52 | 39.29 | 51.76 | Good for reef-level/seasonal interpretation |
| Aquavila | 2022 | Dry | 20 | 35.95 | 29.58 | 42.32 | Good for reef-level/seasonal interpretation |
| Aquavila | 2022 | Rainy | 20 | 49.45 | 40.18 | 58.72 | Good for reef-level/seasonal interpretation |
| Aquavila | 2023 | Dry | 20 | 42.85 | 33.79 | 51.91 | Good for reef-level/seasonal interpretation |
| Aquavila | 2023 | Rainy | 20 | 43.70 | 32.02 | 55.38 | Good for reef-level/seasonal interpretation |
| Aquavila | 2024 | Dry | 20 | 35.25 | 25.36 | 45.14 | Good for reef-level/seasonal interpretation |
| Aquavila | 2024 | Rainy | 20 | 31.80 | 24.58 | 39.02 | Good for reef-level/seasonal interpretation |
| Aquavila | 2025 | Dry | 30 | 33.17 | 25.07 | 41.26 | Good for reef-level/seasonal interpretation |
| Aquavila | 2025 | Rainy | 30 | 42.23 | 31.09 | 53.37 | Good for reef-level/seasonal interpretation |
| Áries | 2021 | Dry | 16 | 10.25 | 1.60 | 18.90 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2022 | Dry | 12 | 0.00 | 0.00 | 0.00 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2023 | Dry | 2 | 0.00 | 0.00 | 0.00 | Low sample size; avoid strong interpretation |
| Áries | 2023 | Rainy | 2 | 0.00 | 0.00 | 0.00 | Low sample size; avoid strong interpretation |
| Áries | 2024 | Dry | 10 | 5.00 | 0.00 | 10.52 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2024 | Rainy | 10 | 19.30 | 6.23 | 32.37 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2025 | Dry | 10 | 26.30 | 4.66 | 47.94 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2025 | Rainy | 10 | 45.90 | 19.13 | 72.67 | Usable descriptively; interpret seasonality cautiously |
| Goiabal | 2021 | Dry | 4 | 16.75 | 0.00 | 39.57 | Low sample size; avoid strong interpretation |
| Goiabal | 2022 | Dry | 2 | 38.00 | 0.00 | 126.94 | Low sample size; avoid strong interpretation |
| Goiabal | 2023 | Dry | 2 | 0.00 | 0.00 | 0.00 | Low sample size; avoid strong interpretation |
| Goiabal | 2023 | Rainy | 2 | 45.50 | 39.15 | 51.85 | Low sample size; avoid strong interpretation |
| Goiabal | 2024 | Dry | 5 | 25.20 | 4.05 | 46.35 | Low sample size; avoid strong interpretation |
| Goiabal | 2024 | Rainy | 5 | 20.80 | 0.00 | 61.97 | Low sample size; avoid strong interpretation |
| Goiabal | 2025 | Dry | 5 | 39.60 | 0.00 | 90.87 | Low sample size; avoid strong interpretation |
| Goiabal | 2025 | Rainy | 5 | 55.20 | 0.00 | 145.68 | Low sample size; avoid strong interpretation |
| Jacarequara | 2021 | Dry | 40 | 38.33 | 32.47 | 44.18 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2022 | Dry | 40 | 15.53 | 12.18 | 18.87 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2023 | Dry | 20 | 32.35 | 25.48 | 39.22 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2023 | Rainy | 20 | 19.90 | 14.34 | 25.46 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2024 | Dry | 20 | 17.35 | 11.72 | 22.98 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2024 | Rainy | 20 | 28.00 | 17.46 | 38.54 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2025 | Dry | 55 | 30.29 | 25.97 | 34.61 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2025 | Rainy | 55 | 62.40 | 52.66 | 72.14 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2021 | Dry | 40 | 36.77 | 28.99 | 44.56 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2022 | Dry | 20 | 14.95 | 10.39 | 19.51 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2022 | Rainy | 20 | 18.40 | 11.82 | 24.98 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2023 | Dry | 20 | 11.90 | 6.56 | 17.24 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2023 | Rainy | 20 | 32.30 | 24.86 | 39.74 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2024 | Dry | 20 | 27.10 | 16.43 | 37.77 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2024 | Rainy | 20 | 38.00 | 27.63 | 48.37 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2025 | Dry | 50 | 24.50 | 14.98 | 34.02 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2025 | Rainy | 50 | 27.06 | 18.65 | 35.47 | Good for reef-level/seasonal interpretation |
| Marauá | 2021 | Dry | 40 | 52.45 | 42.28 | 62.62 | Good for reef-level/seasonal interpretation |
| Marauá | 2022 | Dry | 20 | 26.35 | 21.86 | 30.84 | Good for reef-level/seasonal interpretation |
| Marauá | 2022 | Rainy | 20 | 32.15 | 25.16 | 39.14 | Good for reef-level/seasonal interpretation |
| Marauá | 2023 | Dry | 20 | 26.35 | 22.81 | 29.89 | Good for reef-level/seasonal interpretation |
| Marauá | 2023 | Rainy | 20 | 29.55 | 24.24 | 34.86 | Good for reef-level/seasonal interpretation |
| Marauá | 2024 | Dry | 20 | 44.50 | 38.32 | 50.68 | Good for reef-level/seasonal interpretation |
| Marauá | 2024 | Rainy | 20 | 45.00 | 36.44 | 53.56 | Good for reef-level/seasonal interpretation |
| Marauá | 2025 | Dry | 35 | 46.69 | 38.56 | 54.81 | Good for reef-level/seasonal interpretation |
| Marauá | 2025 | Rainy | 45 | 68.16 | 54.53 | 81.78 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2021 | Dry | 40 | 53.12 | 39.55 | 66.70 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2022 | Dry | 40 | 34.67 | 29.41 | 39.94 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2023 | Dry | 20 | 41.35 | 33.13 | 49.57 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2023 | Rainy | 20 | 47.35 | 39.33 | 55.37 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2024 | Dry | 20 | 54.00 | 41.89 | 66.11 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2024 | Rainy | 20 | 51.75 | 40.59 | 62.91 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2025 | Dry | 20 | 83.75 | 53.94 | 113.56 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2025 | Rainy | 20 | 98.20 | 71.03 | 125.37 | Good for reef-level/seasonal interpretation |
| Romana | 2021 | Dry | 38 | 33.58 | 24.56 | 42.60 | Good for reef-level/seasonal interpretation |
| Romana | 2022 | Dry | 40 | 16.75 | 12.41 | 21.09 | Good for reef-level/seasonal interpretation |
| Romana | 2023 | Dry | 20 | 13.30 | 7.85 | 18.75 | Good for reef-level/seasonal interpretation |
| Romana | 2023 | Rainy | 20 | 22.05 | 13.81 | 30.29 | Good for reef-level/seasonal interpretation |
| Romana | 2024 | Dry | 20 | 30.05 | 18.69 | 41.41 | Good for reef-level/seasonal interpretation |
| Romana | 2024 | Rainy | 20 | 27.55 | 19.64 | 35.46 | Good for reef-level/seasonal interpretation |
| Romana | 2025 | Dry | 35 | 44.89 | 34.56 | 55.21 | Good for reef-level/seasonal interpretation |
| Romana | 2025 | Rainy | 35 | 41.94 | 32.03 | 51.86 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2021 | Dry | 50 | 42.70 | 33.54 | 51.86 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2022 | Dry | 50 | 24.96 | 20.47 | 29.45 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2023 | Dry | 25 | 28.60 | 24.03 | 33.17 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2023 | Rainy | 25 | 25.08 | 20.25 | 29.91 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2024 | Dry | 25 | 24.40 | 17.66 | 31.14 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2024 | Rainy | 25 | 13.96 | 9.34 | 18.58 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2025 | Dry | 75 | 17.39 | 13.14 | 21.63 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2025 | Rainy | 120 | 33.42 | 28.01 | 38.84 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2021 | Dry | 11 | 34.18 | 22.10 | 46.26 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2022 | Dry | 5 | 0.00 | 0.00 | 0.00 | Low sample size; avoid strong interpretation |
| Tio Oscar | 2022 | Rainy | 10 | 32.60 | 23.92 | 41.28 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2023 | Dry | 20 | 18.60 | 14.74 | 22.46 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2024 | Dry | 10 | 22.10 | 13.21 | 30.99 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2024 | Rainy | 10 | 33.10 | 24.67 | 41.53 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2025 | Dry | 20 | 20.70 | 12.28 | 29.12 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2025 | Rainy | 20 | 50.80 | 35.41 | 66.19 | Good for reef-level/seasonal interpretation |
2. Mean oyster density by reef and year
| Reef | Year | Quadrats | Mean density | Lower 95% CI | Upper 95% CI | Mean live density | Mean dead density | Interpretation |
|---|---|---|---|---|---|---|---|---|
| Água Boa | 2021 | 16 | 33.69 | 21.80 | 45.58 | 20.31 | 13.38 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2022 | 30 | 16.97 | 10.87 | 23.06 | 13.13 | 3.83 | Good for reef-level/seasonal interpretation |
| Água Boa | 2023 | 30 | 11.80 | 7.78 | 15.82 | 7.93 | 3.87 | Good for reef-level/seasonal interpretation |
| Água Boa | 2024 | 30 | 13.73 | 8.91 | 18.55 | 3.10 | 10.63 | Good for reef-level/seasonal interpretation |
| Água Boa | 2025 | 50 | 7.08 | 3.56 | 10.60 | 4.72 | 2.36 | Good for reef-level/seasonal interpretation |
| Aquavila | 2021 | 40 | 45.52 | 39.29 | 51.76 | 34.02 | 11.50 | Good for reef-level/seasonal interpretation |
| Aquavila | 2022 | 40 | 42.70 | 36.91 | 48.49 | 33.58 | 9.12 | Good for reef-level/seasonal interpretation |
| Aquavila | 2023 | 40 | 43.27 | 36.22 | 50.33 | 33.80 | 9.47 | Good for reef-level/seasonal interpretation |
| Aquavila | 2024 | 40 | 33.52 | 27.66 | 39.39 | 28.50 | 5.03 | Good for reef-level/seasonal interpretation |
| Aquavila | 2025 | 60 | 37.70 | 30.92 | 44.48 | 29.90 | 7.80 | Good for reef-level/seasonal interpretation |
| Áries | 2021 | 16 | 10.25 | 1.60 | 18.90 | 4.44 | 5.81 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2022 | 12 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2023 | 4 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | Low sample size; avoid strong interpretation |
| Áries | 2024 | 20 | 12.15 | 4.90 | 19.40 | 10.45 | 1.70 | Good for reef-level/seasonal interpretation |
| Áries | 2025 | 20 | 36.10 | 19.90 | 52.30 | 19.45 | 16.65 | Good for reef-level/seasonal interpretation |
| Goiabal | 2021 | 4 | 16.75 | 0.00 | 39.57 | 11.50 | 5.25 | Low sample size; avoid strong interpretation |
| Goiabal | 2022 | 2 | 38.00 | 0.00 | 126.94 | 22.50 | 15.50 | Low sample size; avoid strong interpretation |
| Goiabal | 2023 | 4 | 22.75 | 0.00 | 64.56 | 10.75 | 12.00 | Low sample size; avoid strong interpretation |
| Goiabal | 2024 | 10 | 23.00 | 5.15 | 40.85 | 21.60 | 1.40 | Usable descriptively; interpret seasonality cautiously |
| Goiabal | 2025 | 10 | 47.40 | 7.03 | 87.77 | 31.20 | 16.20 | Usable descriptively; interpret seasonality cautiously |
| Jacarequara | 2021 | 40 | 38.33 | 32.47 | 44.18 | 23.48 | 14.85 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2022 | 40 | 15.53 | 12.18 | 18.87 | 10.55 | 4.97 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2023 | 40 | 26.12 | 21.45 | 30.80 | 16.60 | 9.53 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2024 | 40 | 22.68 | 16.72 | 28.63 | 11.52 | 11.15 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2025 | 110 | 46.35 | 40.28 | 52.41 | 30.39 | 15.95 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2021 | 40 | 36.77 | 28.99 | 44.56 | 26.88 | 9.90 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2022 | 40 | 16.68 | 12.82 | 20.53 | 13.10 | 3.58 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2023 | 40 | 22.10 | 16.63 | 27.57 | 15.48 | 6.62 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2024 | 40 | 32.55 | 25.24 | 39.86 | 23.92 | 8.62 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2025 | 100 | 25.78 | 19.53 | 32.03 | 14.95 | 10.83 | Good for reef-level/seasonal interpretation |
| Marauá | 2021 | 40 | 52.45 | 42.28 | 62.62 | 42.80 | 9.65 | Good for reef-level/seasonal interpretation |
| Marauá | 2022 | 40 | 29.25 | 25.18 | 33.32 | 20.92 | 8.32 | Good for reef-level/seasonal interpretation |
| Marauá | 2023 | 40 | 27.95 | 24.86 | 31.04 | 20.88 | 7.08 | Good for reef-level/seasonal interpretation |
| Marauá | 2024 | 40 | 44.75 | 39.71 | 49.79 | 38.27 | 6.47 | Good for reef-level/seasonal interpretation |
| Marauá | 2025 | 80 | 58.76 | 50.14 | 67.38 | 44.53 | 14.24 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2021 | 40 | 53.12 | 39.55 | 66.70 | 37.58 | 15.55 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2022 | 40 | 34.67 | 29.41 | 39.94 | 24.62 | 10.05 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2023 | 40 | 44.35 | 38.79 | 49.91 | 30.12 | 14.22 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2024 | 40 | 52.88 | 45.01 | 60.74 | 45.60 | 7.28 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2025 | 40 | 90.97 | 71.60 | 110.35 | 61.33 | 29.65 | Good for reef-level/seasonal interpretation |
| Romana | 2021 | 38 | 33.58 | 24.56 | 42.60 | 17.34 | 16.24 | Good for reef-level/seasonal interpretation |
| Romana | 2022 | 40 | 16.75 | 12.41 | 21.09 | 9.43 | 7.32 | Good for reef-level/seasonal interpretation |
| Romana | 2023 | 40 | 17.68 | 12.75 | 22.60 | 13.52 | 4.15 | Good for reef-level/seasonal interpretation |
| Romana | 2024 | 40 | 28.80 | 22.18 | 35.42 | 25.15 | 3.65 | Good for reef-level/seasonal interpretation |
| Romana | 2025 | 70 | 43.41 | 36.43 | 50.40 | 30.49 | 12.93 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2021 | 50 | 42.70 | 33.54 | 51.86 | 28.38 | 14.32 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2022 | 50 | 24.96 | 20.47 | 29.45 | 17.86 | 7.10 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2023 | 50 | 26.84 | 23.60 | 30.08 | 18.22 | 8.62 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2024 | 50 | 19.18 | 14.97 | 23.39 | 10.92 | 8.26 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2025 | 195 | 27.26 | 23.41 | 31.10 | 17.99 | 9.27 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2021 | 11 | 34.18 | 22.10 | 46.26 | 21.82 | 12.36 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2022 | 15 | 21.73 | 11.41 | 32.06 | 15.20 | 6.53 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2023 | 20 | 18.60 | 14.74 | 22.46 | 12.30 | 6.30 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2024 | 20 | 27.60 | 21.49 | 33.71 | 23.30 | 4.30 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2025 | 40 | 35.75 | 26.07 | 45.43 | 24.35 | 11.40 | Good for reef-level/seasonal interpretation |
3. Seasonal density: rainy vs dry periods
The seasonal analysis is useful, but should be treated as complementary because there are some reef-year-season combinations with low sample sizes.
4. Live and dead oyster proportions
Here, proportions are calculated as pooled oyster-level proportions within each reef-year or reef-year-season group. I think this is much easier to explain than the mean of quadrat-level proportions.
| Reef | Year | Quadrats | Live oysters | Dead oysters | Total oysters | Live proportion | Lower 95% CI | Upper 95% CI | Interpretation |
|---|---|---|---|---|---|---|---|---|---|
| Água Boa | 2021 | 16 | 325 | 214 | 539 | 0.603 | 0.561 | 0.643 | Usable descriptively; interpret seasonality cautiously |
| Água Boa | 2022 | 30 | 394 | 115 | 509 | 0.774 | 0.736 | 0.808 | Good for reef-level/seasonal interpretation |
| Água Boa | 2023 | 30 | 238 | 116 | 354 | 0.672 | 0.622 | 0.719 | Good for reef-level/seasonal interpretation |
| Água Boa | 2024 | 30 | 93 | 319 | 412 | 0.226 | 0.188 | 0.269 | Good for reef-level/seasonal interpretation |
| Água Boa | 2025 | 50 | 236 | 118 | 354 | 0.667 | 0.616 | 0.714 | Good for reef-level/seasonal interpretation |
| Aquavila | 2021 | 40 | 1361 | 460 | 1821 | 0.747 | 0.727 | 0.767 | Good for reef-level/seasonal interpretation |
| Aquavila | 2022 | 40 | 1343 | 365 | 1708 | 0.786 | 0.766 | 0.805 | Good for reef-level/seasonal interpretation |
| Aquavila | 2023 | 40 | 1352 | 379 | 1731 | 0.781 | 0.761 | 0.800 | Good for reef-level/seasonal interpretation |
| Aquavila | 2024 | 40 | 1140 | 201 | 1341 | 0.850 | 0.830 | 0.868 | Good for reef-level/seasonal interpretation |
| Aquavila | 2025 | 60 | 1794 | 468 | 2262 | 0.793 | 0.776 | 0.809 | Good for reef-level/seasonal interpretation |
| Áries | 2021 | 16 | 71 | 93 | 164 | 0.433 | 0.359 | 0.509 | Usable descriptively; interpret seasonality cautiously |
| Áries | 2022 | 12 | 0 | 0 | 0 | NA | NA | NA | Usable descriptively; interpret seasonality cautiously |
| Áries | 2023 | 4 | 0 | 0 | 0 | NA | NA | NA | Low sample size; avoid strong interpretation |
| Áries | 2024 | 20 | 209 | 34 | 243 | 0.860 | 0.811 | 0.898 | Good for reef-level/seasonal interpretation |
| Áries | 2025 | 20 | 389 | 333 | 722 | 0.539 | 0.502 | 0.575 | Good for reef-level/seasonal interpretation |
| Goiabal | 2021 | 4 | 46 | 21 | 67 | 0.687 | 0.568 | 0.785 | Low sample size; avoid strong interpretation |
| Goiabal | 2022 | 2 | 45 | 31 | 76 | 0.592 | 0.480 | 0.696 | Low sample size; avoid strong interpretation |
| Goiabal | 2023 | 4 | 43 | 48 | 91 | 0.473 | 0.373 | 0.574 | Low sample size; avoid strong interpretation |
| Goiabal | 2024 | 10 | 216 | 14 | 230 | 0.939 | 0.900 | 0.963 | Usable descriptively; interpret seasonality cautiously |
| Goiabal | 2025 | 10 | 312 | 162 | 474 | 0.658 | 0.614 | 0.700 | Usable descriptively; interpret seasonality cautiously |
| Jacarequara | 2021 | 40 | 939 | 594 | 1533 | 0.613 | 0.588 | 0.637 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2022 | 40 | 422 | 199 | 621 | 0.680 | 0.642 | 0.715 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2023 | 40 | 664 | 381 | 1045 | 0.635 | 0.606 | 0.664 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2024 | 40 | 461 | 446 | 907 | 0.508 | 0.476 | 0.541 | Good for reef-level/seasonal interpretation |
| Jacarequara | 2025 | 110 | 3343 | 1755 | 5098 | 0.656 | 0.643 | 0.669 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2021 | 40 | 1075 | 396 | 1471 | 0.731 | 0.708 | 0.753 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2022 | 40 | 524 | 143 | 667 | 0.786 | 0.753 | 0.815 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2023 | 40 | 619 | 265 | 884 | 0.700 | 0.669 | 0.730 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2024 | 40 | 957 | 345 | 1302 | 0.735 | 0.710 | 0.758 | Good for reef-level/seasonal interpretation |
| Lauro Sodré | 2025 | 100 | 1495 | 1083 | 2578 | 0.580 | 0.561 | 0.599 | Good for reef-level/seasonal interpretation |
| Marauá | 2021 | 40 | 1712 | 386 | 2098 | 0.816 | 0.799 | 0.832 | Good for reef-level/seasonal interpretation |
| Marauá | 2022 | 40 | 837 | 333 | 1170 | 0.715 | 0.689 | 0.741 | Good for reef-level/seasonal interpretation |
| Marauá | 2023 | 40 | 835 | 283 | 1118 | 0.747 | 0.721 | 0.771 | Good for reef-level/seasonal interpretation |
| Marauá | 2024 | 40 | 1531 | 259 | 1790 | 0.855 | 0.838 | 0.871 | Good for reef-level/seasonal interpretation |
| Marauá | 2025 | 80 | 3562 | 1139 | 4701 | 0.758 | 0.745 | 0.770 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2021 | 40 | 1503 | 622 | 2125 | 0.707 | 0.688 | 0.726 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2022 | 40 | 985 | 402 | 1387 | 0.710 | 0.686 | 0.733 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2023 | 40 | 1205 | 569 | 1774 | 0.679 | 0.657 | 0.701 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2024 | 40 | 1824 | 291 | 2115 | 0.862 | 0.847 | 0.876 | Good for reef-level/seasonal interpretation |
| Pinheiro | 2025 | 40 | 2453 | 1186 | 3639 | 0.674 | 0.659 | 0.689 | Good for reef-level/seasonal interpretation |
| Romana | 2021 | 38 | 659 | 617 | 1276 | 0.516 | 0.489 | 0.544 | Good for reef-level/seasonal interpretation |
| Romana | 2022 | 40 | 377 | 293 | 670 | 0.563 | 0.525 | 0.600 | Good for reef-level/seasonal interpretation |
| Romana | 2023 | 40 | 541 | 166 | 707 | 0.765 | 0.733 | 0.795 | Good for reef-level/seasonal interpretation |
| Romana | 2024 | 40 | 1006 | 146 | 1152 | 0.873 | 0.853 | 0.891 | Good for reef-level/seasonal interpretation |
| Romana | 2025 | 70 | 2134 | 905 | 3039 | 0.702 | 0.686 | 0.718 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2021 | 50 | 1419 | 716 | 2135 | 0.665 | 0.644 | 0.684 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2022 | 50 | 893 | 355 | 1248 | 0.716 | 0.690 | 0.740 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2023 | 50 | 911 | 431 | 1342 | 0.679 | 0.653 | 0.703 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2024 | 50 | 546 | 413 | 959 | 0.569 | 0.538 | 0.600 | Good for reef-level/seasonal interpretation |
| Terra Amarela | 2025 | 195 | 3508 | 1807 | 5315 | 0.660 | 0.647 | 0.673 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2021 | 11 | 240 | 136 | 376 | 0.638 | 0.589 | 0.685 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2022 | 15 | 228 | 98 | 326 | 0.699 | 0.648 | 0.747 | Usable descriptively; interpret seasonality cautiously |
| Tio Oscar | 2023 | 20 | 246 | 126 | 372 | 0.661 | 0.612 | 0.708 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2024 | 20 | 466 | 86 | 552 | 0.844 | 0.812 | 0.872 | Good for reef-level/seasonal interpretation |
| Tio Oscar | 2025 | 40 | 974 | 456 | 1430 | 0.681 | 0.657 | 0.705 | Good for reef-level/seasonal interpretation |
5. Associated biodiversity and taxonomic richness
These analyses use the 2021–2024 taxonomic_group_counts
sheet. Maybe this should be framed as observed associated biodiversity
or something like that, and not a full biodiversity inventory, because
richness depends on sampling effort, detection, and the taxonomic
resolution used during monitoring.
| Reef | Year | Biodiversity quadrats | Observed associated taxa richness | Mean taxa richness per quadrat | Mean associated taxa abundance per quadrat | Shannon associated taxa index |
|---|---|---|---|---|---|---|
| Aquavila | 2021 | 40 | 9 | 3.775 | 22.325 | 1.485 |
| Aquavila | 2022 | 40 | 9 | 2.100 | 12.725 | 1.303 |
| Aquavila | 2023 | 40 | 8 | 2.650 | 10.425 | 1.504 |
| Aquavila | 2024 | 40 | 9 | 2.075 | 3.800 | 1.828 |
| Goiabal | 2021 | 4 | 2 | 0.750 | 31.250 | 0.047 |
| Goiabal | 2022 | 2 | 3 | 2.000 | 18.500 | 0.333 |
| Goiabal | 2023 | 4 | 1 | 0.500 | 10.500 | 0.000 |
| Goiabal | 2024 | 10 | 4 | 1.000 | 9.000 | 1.086 |
| Jacarequara | 2021 | 40 | 9 | 2.900 | 24.875 | 1.015 |
| Jacarequara | 2022 | 40 | 9 | 1.825 | 26.775 | 1.016 |
| Jacarequara | 2023 | 40 | 10 | 2.875 | 20.850 | 1.466 |
| Jacarequara | 2024 | 40 | 9 | 1.700 | 9.925 | 1.549 |
| Lauro Sodré | 2021 | 40 | 11 | 3.150 | 33.300 | 1.199 |
| Lauro Sodré | 2022 | 40 | 7 | 1.500 | 7.525 | 1.215 |
| Lauro Sodré | 2023 | 40 | 10 | 2.375 | 21.550 | 0.959 |
| Lauro Sodré | 2024 | 40 | 12 | 1.900 | 7.950 | 1.555 |
| Marauá | 2021 | 40 | 9 | 2.400 | 24.100 | 0.983 |
| Marauá | 2022 | 40 | 8 | 2.100 | 10.525 | 1.124 |
| Marauá | 2023 | 40 | 8 | 2.300 | 19.875 | 0.996 |
| Marauá | 2024 | 40 | 11 | 1.975 | 11.325 | 1.284 |
| Pinheiro | 2021 | 40 | 9 | 2.300 | 48.900 | 0.865 |
| Pinheiro | 2022 | 40 | 7 | 1.300 | 14.925 | 0.850 |
| Pinheiro | 2023 | 40 | 8 | 1.300 | 18.975 | 0.933 |
| Pinheiro | 2024 | 40 | 10 | 1.500 | 20.950 | 1.139 |
| Romana | 2021 | 38 | 7 | 1.895 | 26.368 | 0.870 |
| Romana | 2022 | 40 | 10 | 1.475 | 6.925 | 1.276 |
| Romana | 2023 | 40 | 10 | 2.600 | 24.250 | 1.094 |
| Romana | 2024 | 40 | 13 | 1.700 | 10.275 | 1.561 |
| Terra Amarela | 2021 | 50 | 10 | 2.240 | 36.340 | 1.031 |
| Terra Amarela | 2022 | 50 | 11 | 2.740 | 26.580 | 1.260 |
| Terra Amarela | 2023 | 50 | 10 | 2.240 | 19.780 | 1.370 |
| Terra Amarela | 2024 | 50 | 10 | 1.440 | 9.300 | 1.269 |
| Tio Oscar | 2021 | 11 | 7 | 2.000 | 22.818 | 0.600 |
| Tio Oscar | 2022 | 15 | 5 | 1.600 | 12.533 | 0.682 |
| Tio Oscar | 2023 | 20 | 8 | 1.900 | 9.700 | 1.273 |
| Tio Oscar | 2024 | 20 | 8 | 1.150 | 3.950 | 1.820 |
| Água Boa | 2021 | 16 | 7 | 3.312 | 35.875 | 0.981 |
| Água Boa | 2022 | 30 | 8 | 1.800 | 12.967 | 1.120 |
| Água Boa | 2023 | 30 | 9 | 2.333 | 17.867 | 1.302 |
| Água Boa | 2024 | 30 | 11 | 2.367 | 8.700 | 1.430 |
| Áries | 2021 | 16 | 8 | 2.438 | 15.438 | 1.175 |
| Áries | 2022 | 12 | 1 | 0.833 | 0.833 | 0.000 |
| Áries | 2023 | 4 | 0 | 0.000 | 0.000 | NA |
| Áries | 2024 | 20 | 6 | 1.050 | 11.950 | 1.294 |
6. Exploratory biodiversity-condition relationships
These figures explore whether associated biodiversity appears related to reef condition indicators, especially oyster density, live oyster proportion, and size structure. They should be presented as exploratory relationships rather than causal evidence.
Main takeaway: Overall, these plots suggest that oyster
reefs with higher oyster density, higher live oyster proportion, and
larger oysters also tend to support a richer community of associated
taxa. The pattern is not perfectly consistent across all reefs and
years, but it points in an encouraging direction: reefs that appear
stronger in terms of oyster abundance, condition, and size may also
provide better habitat for other organisms. I also think that these
results should not be interpreted as showing that one factor directly
causes the other. For me, a more likely interpretation is that some
reefs may have better underlying environmental or habitat conditions,
which support both healthier oyster populations and a richer associated
community.
7. Oyster size structure as reef-condition context
These plots could be useful for discussing whether denser or more live-dominated reefs also show different size structures.
