1. Summary

This report presents the results of an exploratory analysis and a Before–After Standardized Effect Size Analysis conducted to evaluate key social measures in Honduras (Trust, Social Cohesion, and Self-Efficacy) using data from the FF HHS. The objective was to find insights into how these social dynamics are evolving over time and to assess whether FF interventions may be driving improvements in these areas.

Overall, the findings reflect mixed results and a lack of clear trends, making it difficult to extract broad, meaningful conclusions about social change across communities. While some communities demonstrated positive shifts, others experienced declines or no significant changes.




2. Objective

Request from the FF leadership:

“Find insights on social measures in Honduras (trust, social cohesion, self-efficacy). Specifically, assess if Fish Forever interventions are driving increases in these measures and if there are links to savings club participation”




3. Exploratory Analaysis

The exploratory analysis is organized into two sections. The first, Available HHS Data, presents tables that summarize the household survey coverage in Honduras, disaggregated by province, municipality, and community across multiple years. The second, Selected HHS Questions, examines the response status and distribution of responses for a set of HHS questions aligned with the core social measures of interest: Trust, Social Cohesion, Self-Efficacy, and Savings Club Participation.

Available HHS Data


#Available Data
merged_hhs_hon <- merged_hhs %>% 
  filter(g1_country == "HND")


table_province <- merged_hhs_hon %>%
  group_by(year, merged_hhs_province) %>%
  summarise(n_surveys = n(), .groups = "drop")

table_municipality <- merged_hhs_hon %>%
  group_by(year, merged_hhs_province, merged_hhs_municipality) %>%
  summarise(n_surveys = n(), .groups = "drop")


table_community <- merged_hhs_hon %>%
  group_by(year, merged_hhs_province, merged_hhs_municipality, merged_hhs_community) %>%
  summarise(n_surveys = n(), .groups = "drop")


#Select relevant columns for analysis
merged_hhs_hon <- merged_hhs_hon %>%
  select(
    year,
    province = merged_hhs_province,
    municipality = merged_hhs_municipality,
    community = merged_hhs_community,
    g8_trust_community,
    g8_trust_community_emergency,
    g8_trust_community_neighbors,
    g8_trust_community_rebuilt,
    g8_trust_local_decision,
    g8_trust_regional_decision,
    g8_my_community_ability,
    g8_fishery_benefit_equal,
    g12_agreement_change_fishing_behavior,
    g12_agreement_individual_behavior,
    g7_savings_club_member,
    g7_emergency_savings_club 
  )


# Check missing data
# Define missing value patterns
missing_vals <- c(NA, 
                  "",
                  #"no_dependence", 
                  #"No dependance",
                  "Not Answered")

# Summarize missing and non-missing values per column
missing_summary <- map_dfr(names(merged_hhs_hon), function(col) {
  data <- merged_hhs_hon[[col]]
  missing_count <- sum(is.na(data) | data %in% missing_vals)
  total_count <- length(data)
  valid_count <- total_count - missing_count

  tibble(
    column = col,
    missing = missing_count,
    valid = valid_count,
    total = total_count,
    missing_pct = round(100 * missing_count / total_count, 1)
  )
})


# Province-level Table
table_province %>%
  arrange(merged_hhs_province, year) %>%
  rename(
    Province = merged_hhs_province,
    `Number of Surveys` = n_surveys,
    Year = year
  ) %>%
  arrange(Province) %>%
  select(Province, `Number of Surveys`, Year) %>%
  kable("html", caption = "Household Surveys by Province and Year") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
Household Surveys by Province and Year
Province Number of Surveys Year
Atlántida 288 2023
Atlántida 288 2024
Atlántida 16 2025
Colón 1160 2021
Colón 275 2023
Colón 1169 2024
Cortés 559 2019
Cortés 632 2024
Islas de la Bahía 227 2019
Islas de la Bahía 275 2024


# Format table for display
table_municipality %>%
  arrange(merged_hhs_province, merged_hhs_municipality, year) %>%
  rename(
    Province = merged_hhs_province,
    Municipality = merged_hhs_municipality,
    `Number of Surveys` = n_surveys,
    Year = year
  ) %>%
  arrange(Province, Municipality) %>%
  select(Province, Municipality, `Number of Surveys`, Year) %>%
  kable("html", caption = "Household Surveys by Municipality and Year") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))
Household Surveys by Municipality and Year
Province Municipality Number of Surveys Year
Atlántida El Porvenir 288 2023
Atlántida La Ceiba 288 2024
Atlántida La Másica 15 2025
Atlántida San Francisco 1 2025
Colón Balfate 3 2021
Colón Balfate 275 2023
Colón Iriona 185 2021
Colón Iriona 187 2024
Colón Limón 102 2021
Colón Limón 104 2024
Colón Santa Fé 284 2021
Colón Santa Fé 291 2024
Colón Santa Rosa de Aguan 283 2021
Colón Santa Rosa de Aguan 292 2024
Colón Trujillo 303 2021
Colón Trujillo 295 2024
Cortés Omoa 275 2019
Cortés Omoa 329 2024
Cortés Puerto Cortés 284 2019
Cortés Puerto Cortés 303 2024
Islas de la Bahía Guanaja 227 2019
Islas de la Bahía Guanaja 148 2024
Islas de la Bahía Jose Santos Guardiola 81 2024
Islas de la Bahía Roatán 19 2024
Islas de la Bahía Utila 27 2024


# Summarize years surveyed per community, including the actual years
community_years_surveyed <- table_community %>%
  group_by(merged_hhs_province, merged_hhs_municipality, merged_hhs_community) %>%
  summarise(
    `Years Surveyed` = paste0(
      n_distinct(year), 
      " (", paste(sort(unique(year)), collapse = ", "), ")"
    ),
    .groups = "drop"
  ) %>%
  arrange(merged_hhs_province, merged_hhs_municipality)

# Create the formatted table
# community_years_surveyed %>%
#   rename(
#     Province = merged_hhs_province,
#     Municipality = merged_hhs_municipality,
#     Community = merged_hhs_community
#   ) %>%
#   kable("html", caption = "Years Surveyed per Community (Count and List of Years)") %>%
#   kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))


# Identify communities surveyed in more than one year
multi_year_communities <- table_community %>%
  group_by(merged_hhs_province, merged_hhs_municipality, merged_hhs_community) %>%
  summarise(n_years = n_distinct(year), .groups = "drop") %>%
  filter(n_years > 1)

# Join back to table_community to get survey counts by year
community_multi_year_surveys <- table_community %>%
  inner_join(multi_year_communities, 
             by = c("merged_hhs_province", "merged_hhs_municipality", "merged_hhs_community")) %>%
  arrange(merged_hhs_province, merged_hhs_municipality, merged_hhs_community, year)

# Prepare data with row numbers
table_data <- community_multi_year_surveys %>%
  rename(
    Province = merged_hhs_province,
    Municipality = merged_hhs_municipality,
    Community = merged_hhs_community,
    Year = year,
    `Number of Surveys` = n_surveys
  ) %>%
  select(Province, Municipality, Community, Year, `Number of Surveys`) %>%
  arrange(Community, Year) %>%
  mutate(row_number = row_number())

# Find last row for each community
break_rows <- table_data %>%
  group_by(Community) %>%
  summarise(last_row = max(row_number), .groups = "drop") %>%
  pull(last_row)

# Create the formatted table
table_out <- table_data %>%
  select(-row_number) %>%
  kbl("html", caption = "Communities Surveyed in Multiple Years with Survey Counts per Year") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

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

# Show final table
table_out
Communities Surveyed in Multiple Years with Survey Counts per Year
Province Municipality Community Year Number of Surveys
Cortés Puerto Cortés Bajamar 2019 120
Cortés Puerto Cortés Bajamar 2024 38
Cortés Omoa Barra Motagua 2019 12
Cortés Omoa Barra Motagua 2024 9
Colón Santa Rosa de Aguan Barra de Aguan 2021 9
Colón Santa Rosa de Aguan Barra de Aguan 2024 11
Colón Trujillo Barranco Blanco 2021 3
Colón Trujillo Barranco Blanco 2024 1
Colón Santa Fé Betulia 2021 16
Colón Santa Fé Betulia 2024 16
Colón Trujillo Capiro 2021 40
Colón Trujillo Capiro 2024 53
Colón Iriona Cocalito 2021 21
Colón Iriona Cocalito 2024 22
Colón Santa Rosa de Aguan Col. Las Lomas 2021 23
Colón Santa Rosa de Aguan Col. Las Lomas 2024 22
Colón Trujillo Cristales 2021 20
Colón Trujillo Cristales 2024 59
Colón Iriona Cusuna 2021 53
Colón Iriona Cusuna 2024 55
Islas de la Bahía Guanaja El Cayo 2019 57
Islas de la Bahía Guanaja El Cayo 2024 112
Islas de la Bahía Guanaja El Pelicano 2019 71
Islas de la Bahía Guanaja El Pelicano 2024 19
Colón Santa Fé Guadalupe 2021 51
Colón Santa Fé Guadalupe 2024 54
Colón Trujillo Jericó 2021 57
Colón Trujillo Jericó 2024 57
Cortés Puerto Cortés Las Sabana 2019 17
Cortés Puerto Cortés Las Sabana 2024 6
Colón Limón Limón 2021 102
Colón Limón Limón 2024 104
Colón Santa Fé Manatí 2021 3
Colón Santa Fé Manatí 2024 3
Islas de la Bahía Guanaja Mangrove Bight 2019 52
Islas de la Bahía Guanaja Mangrove Bight 2024 8
Cortés Omoa Masca 2019 44
Cortés Omoa Masca 2024 35
Cortés Omoa Milla 3 2019 12
Cortés Omoa Milla 3 2024 10
Colón Trujillo Moradel 2021 10
Colón Trujillo Moradel 2024 12
Cortés Omoa Muchilena 2019 24
Cortés Omoa Muchilena 2024 19
Colón Santa Fé Plan Grande 2021 14
Colón Santa Fé Plan Grande 2024 19
Colón Trujillo Puerto Castilla 2021 51
Colón Trujillo Puerto Castilla 2024 44
Colón Santa Fé Quinito 2021 12
Colón Santa Fé Quinito 2024 9
Colón Trujillo Rio Negro 2021 16
Colón Trujillo Rio Negro 2024 22
Colón Santa Fé San Antonio 2021 52
Colón Santa Fé San Antonio 2024 54
Colón Iriona San José de la Punta 2021 27
Colón Iriona San José de la Punta 2024 30
Colón Trujillo San Martín 2021 47
Colón Trujillo San Martín 2024 18
Colón Iriona Sangre Laya 2021 28
Colón Iriona Sangre Laya 2024 26
Colón Santa Fé Santa Fe 2021 136
Colón Santa Fé Santa Fe 2024 136
Colón Trujillo Silin 2021 29
Colón Trujillo Silin 2024 29
Cortés Puerto Cortés Travesía 2019 147
Cortés Puerto Cortés Travesía 2024 8

Note: The communities included in this analysis are those with multiple years of household survey data, as shown in the table above. Only communities with at least two distinct survey years and a minimum of 10 valid responses per year were retained, forming the basis for the Before–After Standardized Effect Size Analysis.

Selected HHS Questions

To assess the social outcomes intended by FF, we selected a subset of HHS questions that align with the core social measures of interest: Trust, Social Cohesion, Self-Efficacy, and Savings Club Participation. This section explores the response status (i.e., answered vs. missing) and the distribution of responses for each selected question. The analysis helps us understand both the availability and variability of data, serving as a foundation for the Before–After Standardized Effect Size Analysis that follow.

Trust

  1. Community Trust: In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries.

  2. Neighbors Trust: Generally speaking, most people in neighboring communities can be trusted.

  3. Local Decision Makers: Local decision-makers/local authorities can be trusted to make decisions that benefit the community over their own interests.

  4. Regional Decision Makers: Regional decision-makers/regional authorities can be trusted to make decisions that benefit the community over their own interests.

Social Cohesion

  1. Community Support: My community will support me with funds and help in the case of emergency.

  2. Post-Emergency Recovery: With the help of community members, the community can be rebuilt after emergencies occur.

  3. Community Ability: My community has the ability to manage my fishery effectively to maximize food and profits.

  4. Fishery Benefit Equality: Do you believe you benefit equally from the fishery as other members of the community?

Self-Efficacy

  1. Change Behavior: I am willing to change my fishing behavior.

  2. Individual Behavior: Through my individual fishing behavior, I can make a meaningful contribution to the sustainability of the fish catch.

Savings Club Participation

  1. Are you or a member of your household a member of a savings club?

  2. Does any household member have access to one of the following emergency funds? (Savings Club)

a. Trust

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community = case_when(
    str_detect(tolower(g8_trust_community), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community), "^disagree$|tidak setuju|b\\. tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community), "^neither$|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))



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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: c) Community Trust *",
    caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )



# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community_neighbors))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community_neighbors)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community_neighbors = case_when(
    str_detect(tolower(g8_trust_community_neighbors), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community_neighbors), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community_neighbors), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community_neighbors), "^agree$|d\\. setuju|setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community_neighbors), "^strongly agree$|strongly_agree|e\\. sangat setuju|sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community_neighbors), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))



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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community_neighbors"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: d) Neighbors Trust *",
    caption = "* d) Generally speaking, most people in neighboring communities can be trusted",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_local_decision))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_local_decision)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_local_decision = case_when(
    str_detect(tolower(g8_trust_local_decision), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_local_decision), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_local_decision), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_local_decision), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_local_decision), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_local_decision), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_local_decision"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: a) Local Decision Makers *",
    caption = "* a) Local decision-makers/ local authorities can be trusted to make decisions that benefit the community over their own interests",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_regional_decision))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_regional_decision)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_regional_decision = case_when(
    str_detect(tolower(g8_trust_regional_decision), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_regional_decision), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_regional_decision), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_regional_decision), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_regional_decision), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_regional_decision), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_regional_decision"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: b) Regional Decision Makers *",
    caption = "* b) Regional decision-makers/ regional authorities can be trusted to make decisions that benefit the community over their own interests",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

b. Social Cohesion

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community_emergency))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community_emergency)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community_emergency = case_when(
    str_detect(tolower(g8_trust_community_emergency), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community_emergency), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community_emergency), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community_emergency), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community_emergency), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community_emergency), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community_emergency"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(
    response_status == "Answered",
    response %in% response_levels,
    province %in% c("Atlántida", "Colón")
  ) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(c("Atlántida", "Colón")))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution for Atlántida and Colón",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: e) Community Support *",
    caption = "* e) My community will support me with funds and help in the case of emergency",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_trust_community_rebuilt))
# 
# # Unique values
# unique(merged_hhs_hon$g8_trust_community_rebuilt)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_trust_community_rebuilt = case_when(
    str_detect(tolower(g8_trust_community_rebuilt), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_trust_community_rebuilt), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_trust_community_rebuilt), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_trust_community_rebuilt), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_trust_community_rebuilt), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_trust_community_rebuilt), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_trust_community_rebuilt"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(
    response_status == "Answered",
    response %in% response_levels,
    province %in% c("Atlántida", "Colón")
  ) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(c("Atlántida", "Colón")))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution for Atlántida and Colón",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: f) Post-Emergency Recovery *",
    caption = "* f) With the help of community members, the community can be rebuilt after emergencies occur",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )


# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_my_community_ability))
# 
# # Unique values
# unique(merged_hhs_hon$g8_my_community_ability)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_my_community_ability = case_when(
    str_detect(tolower(g8_my_community_ability), "^strongly disagree$|strongly_disagree|a\\. sangat tidak setuju|sangat tidak setuju") ~ "Strongly disagree",
    str_detect(tolower(g8_my_community_ability), "^disagree$|b\\. tidak setuju|tidak setuju") ~ "Disagree",
    str_detect(tolower(g8_my_community_ability), "^neither$|neither agree nor disagree|netral|c\\. netral") ~ "Neither agree nor disagree",
    str_detect(tolower(g8_my_community_ability), "^agree$|setuju|d\\. setuju") ~ "Agree",
    str_detect(tolower(g8_my_community_ability), "^strongly agree$|strongly_agree|sangat setuju|e\\. sangat setuju") ~ "Strongly agree",
    str_detect(tolower(g8_my_community_ability), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))




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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g8_my_community_ability"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "How much do you agree with: “My community has the ability to manage my fishery effectively to maximize food and profits”?",
    #caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g8_fishery_benefit_equal))
# 
# # Unique values
# unique(merged_hhs_hon$g8_fishery_benefit_equal)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g8_fishery_benefit_equal = case_when(
    str_detect(tolower(g8_fishery_benefit_equal), "^yes$|a\\. ya|^ya$") ~ "Yes",
    str_detect(tolower(g8_fishery_benefit_equal), "^no$|b\\. tidak|^tidak$") ~ "No",
    str_detect(tolower(g8_fishery_benefit_equal), "no_dependence") ~ "I don’t depend on or benefit from the fishery",
    str_detect(tolower(g8_fishery_benefit_equal), "not answered|^na$|tidak tahu|c\\. tidak tahu") ~ NA_character_,
    TRUE ~ NA_character_
  ))





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

# Define missing values explicitly
missing_vals <- c(NA, 
                  "",
                  #"no_dependence", 
                  #"No dependance",
                  "Not Answered")

# Column to analyze
col <- "g8_fishery_benefit_equal"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Yes", "No", "I don’t depend on or benefit from the fishery"
)

color_mapping <- c(
"Yes" = "#6baed6",
    "No" = "#fb6a4a",
    "I don’t depend on or benefit from the fishery" = "#d9d9d9"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Do you believe you benefit equally from the fishery as other members of the community?",
    #caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

c. Self Efficacy

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g12_agreement_change_fishing_behavior))
# 
# # Unique values
# unique(merged_hhs_hon$g12_agreement_change_fishing_behavior)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g12_agreement_change_fishing_behavior = case_when(
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^strongly disagree$|^strongly_disagree$") ~ "Strongly disagree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^disagree$") ~ "Disagree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^neither agree nor disagree$|^neither$") ~ "Neither agree nor disagree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^agree$") ~ "Agree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^strongly agree$|^strongly_agree$") ~ "Strongly agree",
    str_detect(tolower(g12_agreement_change_fishing_behavior), "^$|not answered|^na$") |
      is.na(g12_agreement_change_fishing_behavior) ~ NA_character_,
    TRUE ~ NA_character_
  ))





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

# Define missing values explicitly
missing_vals <- c(NA, "", "Not Answered", "no_dependence", "No dependance")

# Column to analyze
col <- "g12_agreement_change_fishing_behavior"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: i) Change Behavior *",
    caption = "* i) I am willing to change my fishing behavior",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g12_agreement_individual_behavior))
# 
# # Unique values
# unique(merged_hhs_hon$g12_agreement_individual_behavior)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g12_agreement_individual_behavior = case_when(
    str_detect(tolower(g12_agreement_individual_behavior), "^strongly disagree$|^strongly_disagree$") ~ "Strongly disagree",
    str_detect(tolower(g12_agreement_individual_behavior), "^disagree$") ~ "Disagree",
    str_detect(tolower(g12_agreement_individual_behavior), "^neither agree nor disagree$|^neither$") ~ "Neither agree nor disagree",
    str_detect(tolower(g12_agreement_individual_behavior), "^agree$") ~ "Agree",
    str_detect(tolower(g12_agreement_individual_behavior), "^strongly agree$|^strongly_agree$") ~ "Strongly agree",
    str_detect(tolower(g12_agreement_individual_behavior), "^$|not answered|^na$") |
      is.na(g12_agreement_individual_behavior) ~ NA_character_,
    TRUE ~ NA_character_
  ))





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


# Column to analyze
col <- "g12_agreement_individual_behavior"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Strongly agree", "Agree", "Neither agree nor disagree",
  "Disagree", "Strongly disagree"
)

color_mapping <- c(
  "Strongly agree" = "#6baed6",
  "Agree" = "#9ecae1",
  "Neither agree nor disagree" = "#d9d9d9",
  "Disagree" = "#fcae91",
  "Strongly disagree" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Please rate your agreement with the following statement: f) Individual Behavior *",
    caption = "* f) Through my individual fishing behavior, I can make a meaningful contribution to the sustainability of the fish catch",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

d. Savings Club

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g7_savings_club_member))
# 
# # Unique values
# unique(merged_hhs_hon$g7_savings_club_member)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g7_savings_club_member = case_when(
    str_detect(tolower(g7_savings_club_member), "^yes$") ~ "Yes",
    str_detect(tolower(g7_savings_club_member), "^no$") ~ "No",
    TRUE ~ NA_character_
  ))






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


# Column to analyze
col <- "g7_savings_club_member"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Yes", "No"
)

color_mapping <- c(
"Yes" = "#6baed6",
    "No" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(
    response_status == "Answered",
    response %in% response_levels,
    province %in% c("Atlántida", "Colón")
  ) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(c("Atlántida", "Colón")))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution for Atlántida and Colón",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Are you or a member of your household a member of a savings club?",
    # caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

# # Count number of NA values 
# sum(is.na(merged_hhs_hon$g7_emergency_savings_club))
# 
# # Unique values
# unique(merged_hhs_hon$g7_emergency_savings_club)

# Standardize entries
merged_hhs_hon <- merged_hhs_hon %>%
  mutate(g7_emergency_savings_club = case_when(
    str_detect(tolower(g7_emergency_savings_club), "^yes$") ~ "Yes",
    str_detect(tolower(g7_emergency_savings_club), "^no$") ~ "No",
    str_detect(tolower(g7_emergency_savings_club), "not answered") ~ NA_character_,
    TRUE ~ NA_character_
  ))





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


# Column to analyze
col <- "g7_emergency_savings_club"

# Recode data to identify missing/answered
df_example <- merged_hhs_hon %>%
  mutate(
    response_status = if_else(.data[[col]] %in% missing_vals, "Missing", "Answered"),
    response = .data[[col]]
  )

# 1. Horizontal bar plot of Answered vs Missing per province
plot_status_province <- df_example %>%
  group_by(province, response_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(y = fct_rev(province), x = n, fill = response_status)) +
  geom_col(position = "stack") +
  geom_text(aes(label = paste0(n, " (", percent(prop, accuracy = 1), ")")),
            position = position_stack(vjust = 0.5), size = 3.5) +
  scale_fill_manual(values = c("Answered" = "#1f78b4", "Missing" = "#d9d9d9")) +
  labs(
    title = "Response Status by Province",
    y = "Province",
    x = "Count",
    fill = NULL
  ) +
  theme_minimal()

# 2. Response distribution plot by province
response_levels <- c(
  "Yes", "No"
)

color_mapping <- c(
"Yes" = "#6baed6",
    "No" = "#fb6a4a"
)

plot_responses <- df_example %>%
  filter(response_status == "Answered", response %in% response_levels) %>%
  group_by(province, response) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(province) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  mutate(
    response = factor(response, levels = rev(response_levels)),
    province = factor(province, levels = rev(sort(unique(province))))
  ) %>%
  ggplot(aes(x = province, y = prop, fill = response)) +
  geom_col(position = "fill", width = 0.7) +
  geom_text(aes(label = percent(prop, accuracy = 1)),
            position = position_fill(vjust = 0.5),
            color = "black", size = 3) +
  scale_fill_manual(values = color_mapping) +
  coord_flip() +
  labs(
    title = "Response Distribution by Province",
    x = "Province",
    y = "Proportion",
    fill = "Response"
  ) +
  theme_minimal()

# Combine plots with patchwork
(plot_status_province / plot_responses) +
  plot_layout(heights = c(3, 4)) +
  plot_annotation(
    title = "Does any household member have access to one of the following emergency funds? 6. Savings Club",
    #caption = "* c) In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries",
    theme = theme(
      plot.title = element_text(face = "bold.italic", hjust = 0.5, margin = margin(b = 20)),
      plot.caption = element_text(size = 9, hjust = 0)
    )
  )

Note: One part of the analysis objective (assessing potential links between FF social outcomes and savings club participation) could not be addressed. The plots above indicate that savings club participation was extremely limited across the surveyed sample. This low participation restricts our ability to statistically explore relationships between participation and changes in social measures, and limits the strength of any inferences related to this specific program component.




4. Before-After Standardized Effect Size Analysis

Methodology

A Before–After Standardized Effect Size Analysis was selected to assess changes in key social measures (Trust, Social Cohesion, and Self-Efficacy) across multiple coastal communities in Honduras, using data from the FF HHS. This approach was chosen to explore whether these measures have changed over time in communities where FF interventions have been implemented.

To quantify change, we applied a before–after regression method that compares average responses from the earliest and most recent survey years available for each site. Survey responses were standardized by converting ordinal and binary variables to a continuous 0–100 scale, enabling consistent comparison across indicators.

The analysis was limited to communities with at least two distinct survey years and a minimum of 10 valid (non-missing) responses per year for each question of interest. After applying these criteria, the final dataset included only communities with HHS data collected in either 2019 or 2021 (as baseline) and 2024 (as follow-up), resulting in before–after time spans of 3 or 5 years, depending on the site.

We used robust linear regression models with heteroscedasticity-consistent standard errors to estimate standardized effect sizes (β₁) for each question. These effect sizes represent the average change in responses between the “before” and “after” periods. Each estimate is accompanied by a 95% confidence interval and a significance test, allowing identification of statistically significant positive or negative changes. Results are visualized by question and by community to support comparison across sites and social dimensions.


Interpretation of Effect Sizes (β₁)

  • For Ordinal Questions: The effect size (β₁) represents the average change in responses on a 0-100 scale between the “before” and “after” periods. A positive effect size indicates an increase in the perception or agreement level, while a negative effect size suggests a decline.

  • For Binary Questions: The effect size (β₁) reflects the change in the probability of a “Yes” response in the after period. Positive values imply an increase in the probability of a “Yes” response, while negative values indicate a decline.


Limitations

A key limitation of this analysis stems from the underlying data. Specifically, the absence of a control group (i.e., communities not exposed to FF interventions). Without a counterfactual, it is not possible to determine whether observed changes in social measures are attributable to the program or to other external factors such as economic shifts, governance changes, climate variability, or broader social dynamics. As a result, the findings should be interpreted as descriptive of change over time, rather than as evidence of direct program impact.


Results

a. Trust


# Trust question selection
trust_var <- "g8_trust_community"  

# Define Likert-to-100 mapping
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare data
trust_df <- merged_hhs_hon %>%
  select(community, year, trust = all_of(trust_var)) %>%
  filter(trust %in% names(likert_map)) %>%
  mutate(trust_score = recode(trust, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10) %>%  # Drop year-community combos with fewer than 10
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) > 1) %>%  # Drop communities with only 1 year
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression per community
trust_model <- trust_df %>%
  group_by(community) %>%
  do({
    model <- lm(trust_score ~ period, data = .)
    robust_se <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
    tidy(robust_se) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    ),
    fill_color = case_when(
      significance == "Positive (p < 0.05)" ~ "#6baed6",
      significance == "Negative (p < 0.05)" ~ "#fb6a4a",
      TRUE ~ "white"
    )
  )

# Plot
tru_1 <- ggplot(trust_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Community Trust"),
    fill = "Change"
  ) +
  theme_minimal()

# Trust question selection
trust_var <- "g8_trust_community_neighbors"

# Define Likert-to-100 mapping
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare data
trust_df <- merged_hhs_hon %>%
  select(community, year, trust = all_of(trust_var)) %>%
  filter(trust %in% names(likert_map)) %>%
  mutate(trust_score = recode(trust, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10) %>%  # Drop year-community combos with fewer than 10
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) > 1) %>%  # Drop communities with only 1 year
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression per community
trust_model <- trust_df %>%
  group_by(community) %>%
  do({
    model <- lm(trust_score ~ period, data = .)
    robust_se <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
    tidy(robust_se) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    ),
    fill_color = case_when(
      significance == "Positive (p < 0.05)" ~ "#6baed6",
      significance == "Negative (p < 0.05)" ~ "#fb6a4a",
      TRUE ~ "white"
    )
  )

# Plot
tru_2 <- ggplot(trust_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Neighbors Trust"),
    fill = "Change"
  ) +
  theme_minimal()


# Trust question selection
trust_var <- "g8_trust_local_decision"

# Define Likert-to-100 mapping
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare data
trust_df <- merged_hhs_hon %>%
  select(community, year, trust = all_of(trust_var)) %>%
  filter(trust %in% names(likert_map)) %>%
  mutate(trust_score = recode(trust, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10) %>%  # Drop year-community combos with fewer than 10
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) > 1) %>%  # Drop communities with only 1 year
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression per community
trust_model <- trust_df %>%
  group_by(community) %>%
  do({
    model <- lm(trust_score ~ period, data = .)
    robust_se <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
    tidy(robust_se) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    ),
    fill_color = case_when(
      significance == "Positive (p < 0.05)" ~ "#6baed6",
      significance == "Negative (p < 0.05)" ~ "#fb6a4a",
      TRUE ~ "white"
    )
  )

# Plot
tru_3 <- ggplot(trust_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Local Decision Makers"),
    fill = "Change"
  ) +
  theme_minimal()


# Trust question selection
trust_var <- "g8_trust_regional_decision"

# Define Likert-to-100 mapping
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare data
trust_df <- merged_hhs_hon %>%
  select(community, year, trust = all_of(trust_var)) %>%
  filter(trust %in% names(likert_map)) %>%
  mutate(trust_score = recode(trust, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10) %>%  # Drop year-community combos with fewer than 10
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) > 1) %>%  # Drop communities with only 1 year
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression per community
trust_model <- trust_df %>%
  group_by(community) %>%
  do({
    model <- lm(trust_score ~ period, data = .)
    robust_se <- coeftest(model, vcov = vcovHC(model, type = "HC1"))
    tidy(robust_se) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    ),
    fill_color = case_when(
      significance == "Positive (p < 0.05)" ~ "#6baed6",
      significance == "Negative (p < 0.05)" ~ "#fb6a4a",
      TRUE ~ "white"
    )
  )

# Plot
tru_4 <- ggplot(trust_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Regional Decision Makers"),
    fill = "Change"
  ) +
  theme_minimal()

(tru_1 | tru_2) / (tru_3 | tru_4) +
  plot_annotation(
    title = "Before–after standardized effect sizes for Trust by Community",
    theme = theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
    )
  )

The figure presents before–after standardized effect sizes (β₁) for four trust-related questions across multiple communities. Each subplot corresponds to a distinct dimension of trust:

  • Community Trust: In general, most people in my area can be trusted to comply with regulations related to coastal fisheries and fisheries.

  • Neighbors Trust: Generally speaking, most people in neighboring communities can be trusted.

  • Local Decision Makers: Local decision-makers/local authorities can be trusted to make decisions that benefit the community over their own interests.

  • Regional Decision Makers: Regional decision-makers/regional authorities can be trusted to make decisions that benefit the community over their own interests.

Each point represents the estimated effect size for the before–after change in responses to each trust question for a specific community. The horizontal lines indicate 95% confidence intervals. Point color denotes statistical significance.

Key Takeaways

The analysis of trust-related indicators (Community Trust, Neighbors Trust, Local Decision Makers, and Regional Decision Makers) shows mixed results across communities, with no consistent trend emerging across all dimensions.

  • A large number of communities exhibit non-significant changes.

  • Several communities demonstrate significant positive changes in specific dimensions. For example, Moradel, Plan Grande, Jericó, and El Pelícano showed notable gains, especially in trust toward local and regional decision-makers.

  • Conversely, Muchilena, San Martín, and Bajamar stand out for having significant negative changes across most or all trust dimensions.

b. Social Cohesion

# # 1. Choose one of the Social Cohesion questions
# sc_var <- "g8_trust_community_emergency"
# 
# # Likert-style mapping to 0–100
# likert_map <- c(
#   "Strongly disagree" = 0,
#   "Disagree" = 25,
#   "Neither agree nor disagree" = 50,
#   "Agree" = 75,
#   "Strongly agree" = 100
# )
# 
# # Prepare data (excluding NAs, enforcing ≥10 valid per year, ≥2 years)
# sc_df <- merged_hhs_hon %>%
#   select(community, year, sc = all_of(sc_var)) %>%
#   filter(sc %in% names(likert_map)) %>%
#   mutate(sc_score = recode(sc, !!!likert_map)) %>%
#   group_by(community, year) %>%
#   filter(n() >= 10, sum(!is.na(sc_score)) >= 10) %>%
#   ungroup() %>%
#   group_by(community) %>%
#   filter(n_distinct(year) >= 2) %>%
#   mutate(period = case_when(
#     year == min(year) ~ "before",
#     year == max(year) ~ "after"
#   )) %>%
#   filter(!is.na(period)) %>%
#   ungroup() %>%
#   mutate(period = factor(period, levels = c("before", "after")))
# 
# # Run robust regression per community
# sc_model <- sc_df %>%
#   group_by(community) %>%
#   do({
#     model <- lm(sc_score ~ period, data = .)
#     robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
#     tidy(robust) %>%
#       filter(term == "periodafter") %>%
#       mutate(
#         conf.low = estimate - 1.96 * std.error,
#         conf.high = estimate + 1.96 * std.error
#       )
#   }) %>%
#   ungroup() %>%
#   mutate(
#     significance = case_when(
#       p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
#       p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
#       TRUE ~ "Not Significant"
#     )
#   )
# 
# # Plot effect sizes
# ggplot(sc_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
#   geom_point(shape = 21, size = 3, color = "black") +
#   geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
#   scale_fill_manual(values = c(
#     "Positive (p < 0.05)" = "forestgreen",
#     "Negative (p < 0.05)" = "purple",
#     "Not Significant" = "white"
#   )) +
#   coord_flip() +
#   labs(
#     x = "Community", y = "Effect Size (β1)",
#     title = paste("Before–after effect sizes for", sc_var, "by Community"),
#     fill = "Change"
#   ) +
#   theme_minimal()
# 
# # 1. Choose one of the Social Cohesion questions
# sc_var <- "g8_trust_community_emergency"
# 
# # Likert-style mapping to 0–100
# likert_map <- c(
#   "Strongly disagree" = 0,
#   "Disagree" = 25,
#   "Neither agree nor disagree" = 50,
#   "Agree" = 75,
#   "Strongly agree" = 100
# )
# 
# # Prepare data (excluding NAs, enforcing ≥10 valid per year, ≥2 years)
# sc_df <- merged_hhs_hon %>%
#   select(community, year, sc = all_of(sc_var)) %>%
#   filter(sc %in% names(likert_map)) %>%
#   mutate(sc_score = recode(sc, !!!likert_map)) %>%
#   group_by(community, year) %>%
#   filter(n() >= 10, sum(!is.na(sc_score)) >= 10) %>%
#   ungroup() %>%
#   group_by(community) %>%
#   filter(n_distinct(year) >= 2) %>%
#   mutate(period = case_when(
#     year == min(year) ~ "before",
#     year == max(year) ~ "after"
#   )) %>%
#   filter(!is.na(period)) %>%
#   ungroup() %>%
#   mutate(period = factor(period, levels = c("before", "after")))
# 
# # Run robust regression per community
# sc_model <- sc_df %>%
#   group_by(community) %>%
#   do({
#     model <- lm(sc_score ~ period, data = .)
#     robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
#     tidy(robust) %>%
#       filter(term == "periodafter") %>%
#       mutate(
#         conf.low = estimate - 1.96 * std.error,
#         conf.high = estimate + 1.96 * std.error
#       )
#   }) %>%
#   ungroup() %>%
#   mutate(
#     significance = case_when(
#       p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
#       p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
#       TRUE ~ "Not Significant"
#     )
#   )
# 
# # Plot effect sizes
# ggplot(sc_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
#   geom_point(shape = 21, size = 3, color = "black") +
#   geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
#   scale_fill_manual(values = c(
#     "Positive (p < 0.05)" = "forestgreen",
#     "Negative (p < 0.05)" = "purple",
#     "Not Significant" = "white"
#   )) +
#   coord_flip() +
#   labs(
#     x = "Community", y = "Effect Size (β1)",
#     title = paste("Before–after effect sizes for", sc_var, "by Community"),
#     fill = "Change"
#   ) +
#   theme_minimal()




# Could not fulfill:
  # - At least two distinct years of household surveys with
  # - At least 10 non-NA responses to the question in each year
# 1. Choose one of the Social Cohesion questions
sc_var <- "g8_my_community_ability"

# Likert-style mapping to 0–100
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare data (excluding NAs, enforcing ≥10 valid per year, ≥2 years)
sc_df <- merged_hhs_hon %>%
  select(community, year, sc = all_of(sc_var)) %>%
  filter(sc %in% names(likert_map)) %>%
  mutate(sc_score = recode(sc, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10, sum(!is.na(sc_score)) >= 10) %>%
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) >= 2) %>%
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression per community
sc_model <- sc_df %>%
  group_by(community) %>%
  do({
    model <- lm(sc_score ~ period, data = .)
    robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
    tidy(robust) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    )
  )

# Plot effect sizes
soc_1 <- ggplot(sc_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Community Ability"),
    fill = "Change"
  ) +
  theme_minimal()

# Binary mapping
binary_map <- c(
  "Yes" = 100,
  "No" = 0
)

# Prepare binary data
binary_df <- merged_hhs_hon %>%
  select(community, year, sc = g8_fishery_benefit_equal) %>%
  filter(sc %in% names(binary_map)) %>%
  mutate(sc_score = recode(sc, !!!binary_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10, sum(!is.na(sc_score)) >= 10) %>%
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) >= 2) %>%
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run regression
binary_model <- binary_df %>%
  group_by(community) %>%
  do({
    model <- lm(sc_score ~ period, data = .)
    robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
    tidy(robust) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    )
  )

# Plot
soc_2 <- ggplot(binary_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = "Fishery Benefit Equality",
    fill = "Change"
  ) +
  theme_minimal()


(soc_1 | soc_2) +
  plot_annotation(
    title = "Before–after standardized effect sizes for Social Cohesion by Community",
    theme = theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
    )
  )

The figure presents before–after standardized effect sizes (β₁) for two questions related to social cohesion across multiple communities. Each subplot corresponds to a distinct social cohesion:

  • Community Ability: My community has the ability to manage my fishery effectively to maximize food and profits.

  • Fishery Benefit Equality: Do you believe you benefit equally from the fishery as other members of the community?

Questions on Community Support and Post-Emergency Recovery were excluded from the analysis due to insufficient data. Specifically, these questions were not included in one year of the household survey, resulting in fewer than two years with at least 10 responses per community.

Each point represents the estimated effect size for the before–after change in responses to each social cohesion question for a specific community. The horizontal lines indicate 95% confidence intervals. Point color denotes statistical significance.

Key Takeaways

The analysis of social cohesion indicators (Community Ability and Fishery Benefit Equality) reveals mixed results, with Community Ability showing both positive and negative significant changes across communities, while all significant changes in Fishery Benefit Equality were negative.

  • A large number of communities exhibit non-significant changes.

  • Community Ability shows mixed results, with no clear overall trend. Significant positive changes were observed in Col. Las Lomas, Jericó, Rio Negro, Capiro, and Cusuna. In contrast, 12 communities (including Bajamar, Guadalupe, San Antonio, Muchilena, and Cocalito) experienced significant negative changes.

  • Fishery Benefit Equality shows a more consistent pattern: all statistically significant changes were negative. Of the 22 communities analyzed, 9 (such as Bajamar, Betulia, San Martín, Cocalito, and El Cayo) showed significant declines in perceived equality of fishery benefits.

c. Self Efficacy

# Select which self-efficacy variable to analyze
sef_var <- "g12_agreement_change_fishing_behavior"

# Define Likert scale to 0–100 mapping
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare dataset
sef_df <- merged_hhs_hon %>%
  select(community, year, sef = all_of(sef_var)) %>%
  filter(sef %in% names(likert_map)) %>%
  mutate(sef_score = recode(sef, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10, sum(!is.na(sef_score)) >= 10) %>%  # Keep only year-community combos with ≥10 non-NA responses
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) >= 2) %>%  # At least 2 survey years with data
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression by community
sef_model <- sef_df %>%
  group_by(community) %>%
  do({
    model <- lm(sef_score ~ period, data = .)
    robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
    tidy(robust) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    ),
    fill_color = case_when(
      significance == "Positive (p < 0.05)" ~ "#6baed6",
      significance == "Negative (p < 0.05)" ~ "#fb6a4a",
      TRUE ~ "white"
    )
  )

# Plot effect sizes
self_1 <- ggplot(sef_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Change Behavior"),
    fill = "Change"
  ) +
  theme_minimal()

# Select which self-efficacy variable to analyze
sef_var <- "g12_agreement_individual_behavior"

# Define Likert scale to 0–100 mapping
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

# Prepare dataset
sef_df <- merged_hhs_hon %>%
  select(community, year, sef = all_of(sef_var)) %>%
  filter(sef %in% names(likert_map)) %>%
  mutate(sef_score = recode(sef, !!!likert_map)) %>%
  group_by(community, year) %>%
  filter(n() >= 10, sum(!is.na(sef_score)) >= 10) %>%  # Keep only year-community combos with ≥10 non-NA responses
  ungroup() %>%
  group_by(community) %>%
  filter(n_distinct(year) >= 2) %>%  # At least 2 survey years with data
  mutate(period = case_when(
    year == min(year) ~ "before",
    year == max(year) ~ "after"
  )) %>%
  filter(!is.na(period)) %>%
  ungroup() %>%
  mutate(period = factor(period, levels = c("before", "after")))

# Run robust regression by community
sef_model <- sef_df %>%
  group_by(community) %>%
  do({
    model <- lm(sef_score ~ period, data = .)
    robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
    tidy(robust) %>%
      filter(term == "periodafter") %>%
      mutate(
        conf.low = estimate - 1.96 * std.error,
        conf.high = estimate + 1.96 * std.error
      )
  }) %>%
  ungroup() %>%
  mutate(
    significance = case_when(
      p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
      p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
      TRUE ~ "Not Significant"
    ),
    fill_color = case_when(
      significance == "Positive (p < 0.05)" ~ "#6baed6",
      significance == "Negative (p < 0.05)" ~ "#fb6a4a",
      TRUE ~ "white"
    )
  )

# Plot effect sizes
self_2 <- ggplot(sef_model, aes(x = reorder(community, estimate), y = estimate, fill = significance)) +
  geom_point(shape = 21, size = 3, color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  scale_fill_manual(values = c(
    "Positive (p < 0.05)" = "#6baed6",
    "Negative (p < 0.05)" = "#fb6a4a",
    "Not Significant" = "white"
  )) +
  coord_flip() +
  labs(
    x = "Community", y = "Effect Size (β1)",
    title = paste("Individual Behavior"),
    fill = "Change"
  ) +
  theme_minimal()


(self_1 | self_2) +
  plot_annotation(
    title = "Before–after standardized effect sizes for Self Efficacy by Community",
    theme = theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5)
    )
  )

The figure presents before–after standardized effect sizes (β₁) for two questions related to self efficacy across multiple communities. Each subplot corresponds to a distinct dimension of self efficacy:

  • Change Behavior: I am willing to change my fishing behavior.

  • Individual Behavior: Through my individual fishing behavior, I can make a meaningful contribution to the sustainability of the fish catch.

Each point represents the estimated effect size for the before–after change in responses to each self efficacy question for a specific community. The horizontal lines indicate 95% confidence intervals. Point color denotes statistical significance.

Key Takeaways

The analysis of self efficacy indicators (Change Behavior and Individual Behavior) shows mixed results across communities, with no consistent trend emerging across all dimensions.

  • A large number of communities exhibit non-significant changes.

  • Change Behavior shows mixed outcomes. Significant positive effects were observed in Cocalito, Guadalupe, Capiro, and Santa Fe, while Muchilena, El Cayo, Cusuna, Limón, and San José de la Punta experienced significant negative changes.

  • Individual Behavior presents a similarly polarized pattern. Significant positive changes were observed in Cocalito, Capiro, Santa Fe, and Guadalupe, while Muchilena, El Cayo, San José de la Punta, Limón, Bajamar, and Cusuna showed significant declines in perceived individual influence.

d. Community Level

# 1. Define Likert & Binary mappings
likert_map <- c(
  "Strongly disagree" = 0,
  "Disagree" = 25,
  "Neither agree nor disagree" = 50,
  "Agree" = 75,
  "Strongly agree" = 100
)

binary_map <- c(
  "Yes" = 100,
  "No" = 0
)

# 2. Define all question metadata: Label -> (variable name, type)
question_info <- list(
  # Trust
  "Community Trust" = list(var = "g8_trust_community", type = "likert"),
  "Neighbors Trust" = list(var = "g8_trust_community_neighbors", type = "likert"),
  "Local Decision Makers" = list(var = "g8_trust_local_decision", type = "likert"),
  "Regional Decision Makers" = list(var = "g8_trust_regional_decision", type = "likert"),
  # Social Cohesion
  "Community Ability" = list(var = "g8_my_community_ability", type = "likert"),
  "Fishery Benefit Equality" = list(var = "g8_fishery_benefit_equal", type = "binary"),
  # Self-Efficacy
  "Change Behavior" = list(var = "g12_agreement_change_fishing_behavior", type = "likert"),
  "Individual Behavior" = list(var = "g12_agreement_individual_behavior", type = "likert")
)

all_results <- list()

for (q_label in names(question_info)) {
  q_var <- question_info[[q_label]]$var
  q_type <- question_info[[q_label]]$type

  df <- merged_hhs_hon %>%
    select(community, year, response = all_of(q_var)) %>%
    filter(!is.na(response))

  if (q_type == "likert") {
    df <- df %>%
      filter(response %in% names(likert_map)) %>%
      mutate(score = recode(response, !!!likert_map))
  } else if (q_type == "binary") {
    df <- df %>%
      filter(response %in% names(binary_map)) %>%
      mutate(score = recode(response, !!!binary_map))
  }

  df_cleaned <- df %>%
    group_by(community, year) %>%
    filter(n() >= 10, sum(!is.na(score)) >= 10) %>%
    ungroup() %>%
    group_by(community) %>%
    filter(n_distinct(year) >= 2) %>%
    mutate(period = case_when(
      year == min(year) ~ "before",
      year == max(year) ~ "after"
    )) %>%
    filter(!is.na(period)) %>%
    ungroup() %>%
    mutate(period = factor(period, levels = c("before", "after")))

  if (nrow(df_cleaned) == 0) next  # Skip if no valid data

  model_results <- df_cleaned %>%
    group_by(community) %>%
    do({
      model <- lm(score ~ period, data = .)
      robust <- coeftest(model, vcov = vcovHC(model, "HC1"))
      tidy(robust) %>%
        filter(term == "periodafter") %>%
        mutate(
          conf.low = estimate - 1.96 * std.error,
          conf.high = estimate + 1.96 * std.error
        )
    }) %>%
    ungroup() %>%
    mutate(
      question = q_label,
      significance = case_when(
        p.value < 0.05 & estimate > 0 ~ "Positive (p < 0.05)",
        p.value < 0.05 & estimate < 0 ~ "Negative (p < 0.05)",
        TRUE ~ "Not Significant"
      )
    )

  all_results[[q_label]] <- model_results
}

final_df <- bind_rows(all_results)
# 
# unique_communities <- unique(combined_results$community)


# Define question order and group
question_order <- c(
  "Community Trust",
  "Neighbors Trust",
  "Local Decision Makers",
  "Regional Decision Makers",
  "Community Ability",
  "Fishery Benefit Equality",
  "Change Behavior",
  "Individual Behavior"
)

question_group <- data.frame(
  question = question_order,
  group = c(
    rep("Trust", 4),
    rep("Social Cohesion", 2),
    rep("Self-Efficacy", 2)
  )
)

# Loop through communities
unique_communities <- unique(final_df$community)

plots <- list()  # empty list to store plots

for (comm in unique(final_df$community)) {
  comm_data <- final_df %>%
    filter(community == comm) %>%
    mutate(question = factor(question, levels = rev(question_order)))

  present_questions <- levels(droplevels(comm_data$question))
  y_map <- data.frame(
    question = present_questions,
    y = seq_along(present_questions)
  ) %>%
    left_join(question_group, by = "question")

  separator_y <- y_map %>%
    group_by(group) %>%
    summarise(max_y = max(y), .groups = "drop") %>%
    arrange(max_y) %>%
    mutate(y_sep = max_y + 0.5) %>%
    slice(-n())

  p <- ggplot(comm_data, aes(x = estimate, y = question, fill = significance)) +
    geom_hline(data = separator_y, aes(yintercept = y_sep),
               color = "grey50", linetype = "dashed", linewidth = 0.4) +
    geom_point(shape = 21, size = 3, color = "black") +
    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
    scale_fill_manual(values = c(
      "Positive (p < 0.05)" = "#6baed6",
      "Negative (p < 0.05)" = "#fb6a4a",
      "Not Significant" = "white"
    )) +
    labs(
      x = "Effect Size (β1)", y = "Question",
      title = comm,
      fill = "Significance"
    ) +
    theme_minimal(base_size = 8) +
    theme(legend.position = "none")

  plots[[comm]] <- p  # store each plot
}
wrap_plots(plots, ncol = 3)

Each community plot displays the standardized effect sizes (β1) for up to 8 social indicators, grouped into three thematic areas:

  • Trust: Community Trust, Neighbors Trust, Local Decision Makers, Regional Decision Makers

  • Social Cohesion: Community Ability, Fishery Benefit Equality

  • Self-Efficacy: Change Behavior, Individual Behavior

Each point represents the estimated effect size for the before–after change in responses to each social indicator question. The horizontal lines indicate 95% confidence intervals. Point color denotes statistical significance.

Key Takeaways

At the community level, the analysis reveals strong heterogeneity in outcomes. While some communities exhibit consistent and significant positive social change others show negative or stagnant trends.

  • Communities showing broad positive trends: Capiro, Plan Grande, Jericó, and Moradel stand out with multiple significant positive effect sizes.

  • Communities showing consistent negative trends: Bajamar, Cusuna, El Cayo, Muchilena, Puerto Castilla, and San Martín are notable for widespread significant declines across multiple trust, cohesion, and efficacy indicators.

  • Communities with mixed or moderate changes: Most communities show mix and non-significant results across indicators.




5. Conclusion

This analysis of changes in key social measures (Trust, Social Cohesion, and Self-Efficacy) across coastal communities in Honduras reveals no clear or consistent pattern of change across indicators or communities. While a large number of communities showed non-significant changes across most indicators, others exhibited either positive gains or declines.

The results show substantial variation in effect sizes, both across communities and across questions. For example, a question that showed positive change in one community often showed negative or no change in another. Although some changes were statistically significant, they were scattered and inconsistent, making it difficult to draw broad conclusions about program impact or widespread social progress.

At the site level, the analysis highlights strong heterogeneity. A small number of communities (such as Capiro, Plan Grande, Jericó, and Moradel) stood out for showing significant positive changes across multiple indicators, suggesting localized momentum in social outcomes. Conversely, communities such as Bajamar, Muchilena, Cusuna, and El Cayo experienced consistent negative trends, indicating potential challenges in trust, cohesion, or individual empowerment. Most other communities showed a mix of non-significant results.

Overall, the findings reflect mixed results and a lack of clear trends, making it difficult to extract broad, meaningful conclusions about social change across communities.