1. Resilience Factors: Scores & 95% Confidence Intervals

This section summarizes site-level scores for each resilience factor, including the overall mean and 95% confidence intervals. We also show the distribution of scores across sites to illustrate variation and coverage.

The raw resilience factor scores data underlying these results is available here.


# total # of sites (across everything)
total_sites <- resilience_scores %>% distinct(ma_id) %>% nrow()

# CI 95%:
ci_t <- function(x, conf = 0.95) {
  x <- x[is.finite(x)]
  n <- length(x); m <- mean(x); s <- sd(x)
  if (n < 2 || is.na(s) || s == 0) return(c(lcl = NA_real_, ucl = NA_real_))
  tcrit <- qt(1 - (1 - conf)/2, df = n - 1)
  se <- s / sqrt(n)
  c(lcl = m - tcrit * se, ucl = m + tcrit * se)
}

# in case we want to use bootstraping:
ci_boot <- function(x, conf = 0.95, B = 5000) {
  x <- x[is.finite(x)]
  n <- length(x)
  if (n < 2) return(c(lcl = NA_real_, ucl = NA_real_))
  mboot <- replicate(B, mean(sample(x, n, replace = TRUE)))
  qs <- quantile(mboot, probs = c((1-conf)/2, 1 - (1-conf)/2), names = FALSE)
  c(lcl = qs[1], ucl = qs[2])
}

use_bootstrap <- FALSE  # set TRUE if you prefer bootstrap CIs

rf_summary <- resilience_scores %>%
  filter(is.finite(score)) %>%
  group_by(resilience_factor) %>%
  summarise(
    n_sites   = n_distinct(ma_id),                      
    mean      = mean(score),
    lcl_ucl   = list(if (use_bootstrap) ci_boot(score) else ci_t(score)),
    .groups   = "drop"
  ) %>%
  mutate(
    lcl = lcl_ucl %>% map_dbl(1),
    ucl = lcl_ucl %>% map_dbl(2),
    # clamp to 1–4 scale
    lcl = pmax(1, pmin(4, lcl)),
    ucl = pmax(1, pmin(4, ucl)),
    coverage_pct = round(n_sites / total_sites * 100, 0)
  ) %>%
  select(resilience_factor, n_sites, coverage_pct, mean, lcl, ucl) %>%
  mutate(across(c(mean, lcl, ucl, coverage_pct), ~ round(.x, 2))) %>%
  arrange(resilience_factor)

rf_summary
#> # A tibble: 6 × 6
#>   resilience_factor                       n_sites coverage_pct  mean   lcl   ucl
#>   <chr>                                     <int>        <dbl> <dbl> <dbl> <dbl>
#> 1 Adaptive Capacity to Climate Change          12            4  2.58  2.26  2.91
#> 2 Capacity for Collective Action              127           47  3.06  2.93  3.2 
#> 3 Coastal Biodiversity & Ecosystem Servi…     272          100  1.42  1.35  1.5 
#> 4 Coastal Fishery Co-Management               272          100  2.6   2.44  2.75
#> 5 Coastal Fishery Policy and Governance       272          100  3.19  3.15  3.24
#> 6 Sustainable Livelihoods & Food Security     150           55  1.67  1.54  1.81

# set desired order
rf_order <- c(
  "Coastal Fishery Policy and Governance",
  "Coastal Fishery Co-Management",
  "Adaptive Capacity to Climate Change",
  "Capacity for Collective Action",
  "Sustainable Livelihoods & Food Security",
  "Coastal Biodiversity & Ecosystem Services"
)

rf_summary %>%
  mutate(resilience_factor = factor(resilience_factor, levels = rf_order)) %>%
  arrange(resilience_factor) %>%
  select(
    `Resilience Factor` = resilience_factor,
    `Number of Sites`   = n_sites,
    `Coverage (%)`      = coverage_pct,
    `Mean Score`        = mean,
    `95% CI (Lower)`    = lcl,
    `95% CI (Upper)`    = ucl
  ) %>%
  mutate(
    `Coverage (%)`   = round(`Coverage (%)`, 1),
    `Mean Score`     = round(`Mean Score`, 2),
    `95% CI (Lower)` = round(`95% CI (Lower)`, 2),
    `95% CI (Upper)` = round(`95% CI (Upper)`, 2)
  ) %>%
  kbl(
    caption = "Summary of Resilience Factors with Mean Scores and 95% Confidence Intervals",
    align = "lrrrrr",
    booktabs = TRUE
  ) %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  column_spec(1, bold = TRUE) %>%
  add_header_above(c(" " = 3, "Resilience Scores" = 3))
Summary of Resilience Factors with Mean Scores and 95% Confidence Intervals
Resilience Scores
Resilience Factor Number of Sites Coverage (%) Mean Score 95% CI (Lower) 95% CI (Upper)
Coastal Fishery Policy and Governance 272 100 3.19 3.15 3.24
Coastal Fishery Co-Management 272 100 2.60 2.44 2.75
Adaptive Capacity to Climate Change 12 4 2.58 2.26 2.91
Capacity for Collective Action 127 47 3.06 2.93 3.20
Sustainable Livelihoods & Food Security 150 55 1.67 1.54 1.81
Coastal Biodiversity & Ecosystem Services 272 100 1.42 1.35 1.50 


# 1) Desired order (top → bottom after coord_flip)
rf_order <- c(
  "Coastal Fishery Policy and Governance",
  "Coastal Fishery Co-Management",
  "Adaptive Capacity to Climate Change",
  "Capacity for Collective Action",
  "Sustainable Livelihoods & Food Security",
  "Coastal Biodiversity & Ecosystem Services"
)

# 2) One row per site × factor
rf_site_factor <- resilience_scores %>%
  filter(is.finite(score)) %>%
  mutate(resilience_factor = str_squish(resilience_factor)) %>%
  group_by(ma_id, ma_name, resilience_factor) %>%
  summarise(score = mean(score), .groups = "drop") %>%
  mutate(resilience_factor = fct_relevel(resilience_factor, rf_order))

# 3) Summary df must use the SAME levels
rf_summary_plot <- rf_summary %>%
  mutate(resilience_factor = str_squish(resilience_factor)) %>%
  mutate(resilience_factor = fct_relevel(resilience_factor, rf_order))

# 4) Plot (reverse limits because of coord_flip)
ggplot(rf_site_factor, aes(x = resilience_factor, y = score)) +
  geom_violin(
    aes(fill = resilience_factor, color = resilience_factor),
    trim = FALSE, alpha = 0.25, linewidth = 0.7
  ) +
  geom_jitter(
    aes(color = resilience_factor),
    width = 0.25, height = 0.25, alpha = 0.55, size = 1.6, show.legend = FALSE
  ) +
  geom_segment(
    data = rf_summary_plot,
    aes(x = resilience_factor, xend = resilience_factor, y = lcl, yend = ucl),
    color = "black", linewidth = 0.8, inherit.aes = FALSE
  ) +
  geom_point(
    data = rf_summary_plot,
    aes(x = resilience_factor, y = mean),
    color = "black", size = 2.6, inherit.aes = FALSE
  ) +
  geom_text(
    data = rf_summary_plot,
    aes(x = resilience_factor, y = mean, label = paste0("n=", n_sites)),
    vjust = -1.2, size = 3, color = "black", inherit.aes = FALSE
  ) +
  scale_x_discrete(limits = rev(rf_order)) +  # enforce order with coord_flip
  coord_flip() +
  scale_y_continuous(limits = c(1, 4), breaks = 1:4) +
  scale_fill_brewer(palette = "Set2", guide = "none") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  labs(
    title = "Site-level score distribution with means and 95% CIs",
    x = NULL, y = "Score (1–4)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank()
  )



3. Years of Program Engagement by CPUE Status

In this section we combine CPUE trend classifications with the year when the management body was formed. We calculate the number of years of program engagement for each site and report the average years across sites with Increasing, Stable, and Decreasing CPUE. This helps explore whether duration of engagement is linked to CPUE outcomes.


# 1) One mgmt_formed per site (take earliest if duplicates) → parse m/d/Y
rf_mgmt <- resilience_scores %>%
  distinct(ma_id, ma_name, mgmt_formed) %>%
  mutate(mgmt_formed = mdy(mgmt_formed)) %>%
  group_by(ma_id, ma_name) %>%
  summarise(mgmt_start = min(mgmt_formed, na.rm = TRUE), .groups = "drop") %>%
  mutate(mgmt_year = year(mgmt_start))

# 2) Make sure your CPUE table carries ma_id
#    If cpue_latest_vs_hist currently has only `site` (ma_name), re-create with id:
site_year_cpue <- cpue_trip_rf %>%
  group_by(ma_id, ma_name, year) %>%
  summarise(avg_cpue = mean(cpue_kg_per_trip, na.rm = TRUE), .groups = "drop")

cpue_latest_vs_hist <- site_year_cpue %>%
  group_by(ma_id, ma_name) %>%
  summarise(
    last_year        = max(year, na.rm = TRUE),
    latest_average   = mean(avg_cpue[year == max(year, na.rm = TRUE)], na.rm = TRUE),
    historic_average = if (n_distinct(year) > 1)
                         mean(avg_cpue[year < max(year, na.rm = TRUE)], na.rm = TRUE)
                       else NA_real_,
    .groups = "drop"
  ) %>%
  mutate(site = ma_name)

# 3) Classify status (use your tolerance)
tol <- 0.15
cpue_status <- cpue_latest_vs_hist %>%
  mutate(
    status = case_when(
      is.na(historic_average) ~ "only one year of data",
      latest_average >= (1 + tol) * historic_average ~ "increasing",
      latest_average <= (1 - tol) * historic_average ~ "decreasing",
      TRUE ~ "stable"
    )
  )

# 4) Join management start year and compute years worked
cpue_with_mgmt <- cpue_status %>%
  left_join(rf_mgmt, by = c("ma_id", "ma_name")) %>%
  mutate(
    years_worked = if_else(!is.na(mgmt_year), last_year - mgmt_year, NA_integer_),
    years_worked = if_else(years_worked < 0, NA_integer_, years_worked)  # guard bad dates
    # If you want inclusive counting, use: last_year - mgmt_year + 1
  )

# 5) Averages of years worked by status (exclude 1-year sites)
years_worked_by_status <- cpue_with_mgmt %>%
  filter(status %in% c("increasing", "stable", "decreasing")) %>%
  summarise(
    n_sites      = n(),
    mean_years   = mean(years_worked, na.rm = TRUE),
    median_years = median(years_worked, na.rm = TRUE),
    sd_years     = sd(years_worked, na.rm = TRUE),
    .by = status
  ) %>%
  arrange(match(status, c("increasing", "stable", "decreasing")))

years_worked_by_status %>%
  mutate(Status = stringr::str_to_title(status)) %>%  # Capitalize
  select(
    Status,
    `Number of Sites`            = n_sites,
    `Mean Years` = mean_years,
    `Median Years` = median_years,
    `SD (Years)`                 = sd_years
  ) %>%
  mutate(
    `Mean Years`   = round(`Mean Years`, 1),
    `Median Years` = round(`Median Years`, 1),
    `SD (Years)`                   = round(`SD (Years)`, 1)
  ) %>%
  kbl(
    caption = "Years Since Management Body Formation by CPUE Trend Status",
    align = "lrrrr",
    booktabs = TRUE
  ) %>%
  kable_classic(full_width = FALSE, html_font = "Arial") %>%
  column_spec(1, bold = TRUE) %>%
  add_header_above(c(" " = 2, "Years Since Management Body Formation" = 3))
Years Since Management Body Formation by CPUE Trend Status
Years Since Management Body Formation
Status Number of Sites Mean Years Median Years SD (Years)
Increasing 53 8.5 4 8.4
Stable 30 4.5 5 3.0
Decreasing 82 8.8 5 8.4


There is no evident relationship between the number of years a site has been under management and its CPUE trend. Sites with longer engagement histories appear just as likely to experience decreasing CPUE as increasing.