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.
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”
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 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"))
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"))
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
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.
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
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.
Social Cohesion
Community Support: My community will support me with funds and help in the case of emergency.
Post-Emergency Recovery: With the help of community members, the community can be rebuilt after emergencies occur.
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?
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.
Savings Club Participation
Are you or a member of your household a member of a savings club?
Does any household member have access to one of the following emergency funds? (Savings Club)
# # 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)
)
)
# # 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)
)
)
# # 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.
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.
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.
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.
# 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.
# 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.
# 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.
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.