Purpose

This analysis explores household income, income-source dependence, and livelihood resilience in the Mozambique HHS dataset.

The analysis focuses on seven questions:

  1. Explore how reported income sources have changed over time across Mozambique.
  2. Explore whether households with higher income dependence on artisanal fishing or farming report higher or lower total real monthly household income.
  3. Explore dependence on artisanal fishing income, overall and by municipality-year.
  4. Assess whether artisanal fishing income dependence is associated with the two Sustainable Livelihoods indicators: food security and ability to meet household needs.
  5. Repeat the same association analysis for farming income dependence.
  6. Explore whether the two Sustainable Livelihoods conditions are becoming less negative through time across Mozambique municipalities, using food worry and financial strain as diagnostic negative indicators.
  7. Explore CPUE coverage/precision by municipality and test whether municipality-year CPUE is associated with HHS-derived monthly household income from artisanal fishing.

This is an exploratory descriptive analysis. It is not a causal impact evaluation. Patterns should be interpreted as associations in repeated cross-sectional HHS data, with year and site composition changing over time.

Data and methodological choices

Income comparability

The HHS income variable is treated as average monthly household income in nominal MZN. Because all observations are from Mozambique, the main comparison is made in real local currency, not USD.

Income is converted to constant 2021 MZN as:

\[ \text{real monthly household income}_{2021\ MZN} = \text{nominal monthly household income}_{MZN} \times \text{deflator to 2021 MZN} \]

The CPI deflators used are included directly in the code below. The 2025 and 2026 values should be treated as provisional and replaced when final annual CPI values are available.

Income outliers

Income is usually skewed, and a small number of very large values can dominate plots and means. For the main descriptive plots, the analysis removes income observations above the 99th percentile of real income within each survey year.

This keeps the vast majority of records while preventing a few extreme values from controlling the visual scale. Medians are emphasized because they are more robust to skewed income distributions than means.

Income-source percentages

The income-source columns are interpreted as the reported percentage of household income coming from each source. For these income-source variables only, NA is treated as zero when at least one source is reported in the row, because blank cells appear to represent sources not used by the household.

The analysis also checks whether the source percentages sum to approximately 100%. This is important because some early years contain records where reported source percentages sum to less than or more than 100%, which can affect compositional interpretation.

Income dependence and total income

To explore whether income-source dependence is associated with total household income, the analysis compares real monthly household income in constant 2021 MZN across dependence bins for artisanal fishing and farming. It also fits exploratory adjusted log-income models controlling for survey year and municipality. The model coefficients are expressed as the approximate percent difference in total income associated with each 10 percentage point increase in income dependence. These estimates are descriptive associations, not causal effects.

Sustainable Livelihoods indicators

Two HHS-based indicators are used:

  1. Food security: respondent reports they Never worried about not having enough food for everyone in the household during the last 12 months.
  2. Income sufficiency / ability to meet household needs: respondent reports that the household can cover its needs Fairly easy, Easy, or Very easy.
## [1] "/Users/marianoviz/Desktop/R Projects and Stuff/ff_hhs_data_processing/data/raw/all_hhs_moz.csv"
## Rows: 7,297
## Columns: 23
## $ year                                                     <dbl> 2021, 2021, 2…
## $ g8_fishery_benefit_equal                                 <chr> "Yes", "No", …
## $ g8_trust_local_decision                                  <chr> "Neither agre…
## $ g8_my_community_ability                                  <chr> "Agree", "Nei…
## $ g12_agreement_community_participation                    <chr> NA, NA, NA, N…
## $ g11_food_worry                                           <chr> "Sometimes", …
## $ g13_hh_ends_meet                                         <chr> "With great d…
## $ g13_hh_average_income                                    <dbl> 2500, 3000, 5…
## $ g4_hh_average_income_source_a_income_farming             <dbl> 100, 30, 100,…
## $ g4_hh_average_income_source_b_income_harvesting          <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_c_income_fishing_artisanal   <dbl> 0, 0, 0, 0, 6…
## $ g4_hh_average_income_source_d_income_fishing_aquaculture <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_e_income_buying_trading      <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_f_income_processing          <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_g_income_extraction          <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_h_income_tourism             <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_i_income_other_wage          <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_j_income_industrial          <dbl> 0, 0, 0, 0, 0…
## $ g4_hh_average_income_source_k_income_other               <dbl> 0, 70, 0, 0, …
## $ g1_country                                               <chr> "MOZ", "MOZ",…
## $ g1_province                                              <chr> "Inhambane", …
## $ g1_municipality                                          <chr> "Inhassoro", …
## $ g1_community                                             <chr> "Petane1", "P…
# Approximate CI for a proportion.
prop_ci <- function(x) {
  x <- x[!is.na(x)]
  n <- length(x)

  if (n == 0) {
    return(tibble(n = 0, pct = NA_real_, ci_low = NA_real_, ci_high = NA_real_))
  }

  p <- mean(x)
  se <- sqrt(p * (1 - p) / n)

  tibble(
    n = n,
    pct = 100 * p,
    ci_low = 100 * pmax(0, p - 1.96 * se),
    ci_high = 100 * pmin(1, p + 1.96 * se)
  )
}

# Mean CI for a continuous percentage variable.
mean_ci <- function(x) {
  x <- x[!is.na(x)]
  n <- length(x)

  if (n == 0) {
    return(tibble(n = 0, mean = NA_real_, ci_low = NA_real_, ci_high = NA_real_))
  }

  m <- mean(x)
  se <- sd(x) / sqrt(n)

  tibble(
    n = n,
    mean = m,
    ci_low = pmax(0, m - 1.96 * se),
    ci_high = pmin(100, m + 1.96 * se)
  )
}

# Bootstrap CI for a median.
bootstrap_median_ci <- function(x, n_boot = 1000, conf = 0.95) {
  x <- x[!is.na(x)]

  if (length(x) == 0) {
    return(tibble(median = NA_real_, ci_low = NA_real_, ci_high = NA_real_))
  }

  if (length(x) < 5) {
    return(tibble(median = median(x), ci_low = NA_real_, ci_high = NA_real_))
  }

  boot_medians <- replicate(
    n_boot,
    median(sample(x, size = length(x), replace = TRUE), na.rm = TRUE)
  )

  alpha <- 1 - conf

  tibble(
    median = median(x),
    ci_low = quantile(boot_medians, alpha / 2, na.rm = TRUE),
    ci_high = quantile(boot_medians, 1 - alpha / 2, na.rm = TRUE)
  )
}

dependency_bins <- function(x) {
  cut(
    x,
    breaks = c(-0.1, 0, 25, 50, 75, 100),
    labels = c("0%", "1–25%", "26–50%", "51–75%", "76–100%")
  )
}
# CPI deflators to convert nominal MZN to constant 2021 MZN.
cpi_deflators <- tribble(
  ~year, ~cpi_2010_100, ~deflator_to_2021_mzn,
  2019, 182.353, 1.1012,
  2020, 188.706, 1.0641,
  2021, 200.799, 1.0000,
  2022, 221.440, 0.9068,
  2023, 237.222, 0.8465,
  2024, 246.898, 0.8133,
  2025, 257.683, 0.7792,
  2026, 269.021, 0.7464
)

hhs_numeric_sources <- hhs_raw %>%
  mutate(
    year = as.integer(year),
    nominal_monthly_hh_income_mzn = parse_number(as.character(.data[[income_col]]))
  ) %>%
  mutate(
    across(
      all_of(source_cols),
      ~ parse_number(as.character(.x))
    )
  ) %>%
  mutate(
    source_module_available = if_any(all_of(source_cols), ~ !is.na(.x))
  )

# Treat source NAs as zero only when the row has at least one reported source.
hhs_sources_zeroed <- hhs_numeric_sources %>%
  mutate(
    across(
      all_of(source_cols),
      ~ case_when(
        !source_module_available ~ NA_real_,
        is.na(.x) ~ 0,
        TRUE ~ .x
      )
    )
  ) %>%
  mutate(
    across(
      all_of(source_cols),
      ~ pmin(pmax(.x, 0), 100)
    )
  )

hhs <- hhs_sources_zeroed %>%
  left_join(cpi_deflators, by = "year") %>%
  mutate(
    real_monthly_hh_income_2021_mzn =
      nominal_monthly_hh_income_mzn * deflator_to_2021_mzn,

    # Sustainable Livelihoods outcomes.
    food_secure = case_when(
      is.na(g11_food_worry) ~ NA_real_,
      str_to_lower(str_trim(g11_food_worry)) == "never" ~ 1,
      TRUE ~ 0
    ),

    income_sufficient = case_when(
      is.na(g13_hh_ends_meet) ~ NA_real_,
      str_to_lower(str_trim(g13_hh_ends_meet)) %in%
        c("fairly easy", "easy", "very easy") ~ 1,
      TRUE ~ 0
    ),

    food_worry = case_when(
      is.na(g11_food_worry) ~ NA_real_,
      str_to_lower(str_trim(g11_food_worry)) %in% c("sometimes", "often") ~ 1,
      TRUE ~ 0
    ),

    financial_strain = case_when(
      is.na(g13_hh_ends_meet) ~ NA_real_,
      str_to_lower(str_trim(g13_hh_ends_meet)) %in%
        c("with difficulty", "with great difficulty") ~ 1,
      TRUE ~ 0
    ),

    source_total_pct = rowSums(across(all_of(source_cols)), na.rm = TRUE),
    source_total_quality = case_when(
      !source_module_available ~ "No source data",
      source_total_pct >= 95 & source_total_pct <= 105 ~ "Approximately 100%",
      source_total_pct < 95 ~ "Below 95%",
      source_total_pct > 105 ~ "Above 105%"
    ),

    fishing_artisanal_income_pct =
      .data[["g4_hh_average_income_source_c_income_fishing_artisanal"]],

    farming_income_pct =
      .data[["g4_hh_average_income_source_a_income_farming"]]
  )

# Remove income outliers above the 99th percentile within each survey year.
income_cutoffs <- hhs %>%
  filter(
    !is.na(real_monthly_hh_income_2021_mzn),
    real_monthly_hh_income_2021_mzn >= 0
  ) %>%
  group_by(year) %>%
  summarise(
    income_p99_by_year = quantile(real_monthly_hh_income_2021_mzn, 0.99, na.rm = TRUE),
    .groups = "drop"
  )

hhs <- hhs %>%
  left_join(income_cutoffs, by = "year") %>%
  mutate(
    valid_income = !is.na(real_monthly_hh_income_2021_mzn) &
      real_monthly_hh_income_2021_mzn >= 0,

    income_outlier_p99_by_year = valid_income &
      real_monthly_hh_income_2021_mzn > income_p99_by_year,

    real_monthly_hh_income_2021_mzn_clean = if_else(
      valid_income & !income_outlier_p99_by_year,
      real_monthly_hh_income_2021_mzn,
      NA_real_
    )
  )

1. Data coverage and quality checks

HHS coverage by year

coverage_year <- hhs %>%
  count(year, name = "n_hhs") %>%
  arrange(year)

coverage_year %>%
  kable(caption = "HHS records by survey year")
HHS records by survey year
year n_hhs
2019 1460
2021 2493
2023 313
2024 711
2025 1865
2026 455 
ggplot(coverage_year, aes(x = factor(year), y = n_hhs)) +
  geom_col() +
  geom_text(aes(label = n_hhs), vjust = -0.3, size = 3.5) +
  scale_y_continuous(labels = comma, expand = expansion(mult = c(0, 0.08))) +
  labs(
    title = "Mozambique HHS coverage by year",
    x = "Survey year",
    y = "HHS records"
  )

HHS coverage by municipality and year

coverage_municipality_year <- hhs %>%
  count(g1_municipality, year, name = "n_hhs")

ggplot(
  coverage_municipality_year,
  aes(
    x = factor(year),
    y = fct_reorder(g1_municipality, n_hhs, .fun = sum),
    fill = n_hhs
  )
) +
  geom_tile(color = "white") +
  geom_text(aes(label = n_hhs), size = 3) +
  scale_fill_viridis_c(labels = comma, option = "C") +
  labs(
    title = "HHS coverage by municipality and year",
    x = "Survey year",
    y = NULL,
    fill = "HHS records"
  ) +
  theme(panel.grid = element_blank())

Income-source percentage quality check

source_total_quality_year <- hhs %>%
  filter(source_module_available) %>%
  count(year, source_total_quality, name = "n") %>%
  group_by(year) %>%
  mutate(pct = 100 * n / sum(n)) %>%
  ungroup()

source_total_quality_year %>%
  arrange(year, source_total_quality) %>%
  mutate(pct = round(pct, 1)) %>%
  kable(
    caption = "Quality check: do reported income-source percentages sum to approximately 100%?"
  )
Quality check: do reported income-source percentages sum to approximately 100%?
year source_total_quality n pct
2019 Above 105% 89 6.1
2019 Approximately 100% 946 64.8
2019 Below 95% 425 29.1
2021 Above 105% 29 1.2
2021 Approximately 100% 2192 87.9
2021 Below 95% 272 10.9
2023 Approximately 100% 313 100.0
2024 Approximately 100% 711 100.0
2025 Approximately 100% 1865 100.0
2026 Approximately 100% 455 100.0 
# ggplot(
#   hhs %>% filter(source_module_available),
#   aes(x = factor(year), y = source_total_pct)
# ) +
#   geom_hline(yintercept = 100, linetype = "dashed", alpha = 0.6) +
#   geom_boxplot(outlier.alpha = 0.15) +
#   coord_cartesian(ylim = c(0, 150)) +
#   labs(
#     title = "Reported income-source percentages by year",
#     subtitle = "Each row should ideally sum to approximately 100%; the dashed line marks 100%",
#     x = "Survey year",
#     y = "Sum of reported income-source percentages"
#   )

Note on interpretation

The income-source variables appear most internally consistent from 2023 onward, with source shares generally summing to 100%. Some earlier records, especially in 2019 and to a lesser extent 2021, have source percentages that sum below or above 100%. For this reason, the source analysis below emphasizes mean reported share per source rather than treating every row as a perfectly closed composition.

Income outlier removal summary

income_outlier_summary <- hhs %>%
  group_by(year) %>%
  summarise(
    total_records = n(),
    records_with_nonmissing_income = sum(valid_income, na.rm = TRUE),
    income_p99_by_year = first(income_p99_by_year),
    records_removed_as_income_outliers = sum(income_outlier_p99_by_year, na.rm = TRUE),
    pct_removed = 100 * records_removed_as_income_outliers / records_with_nonmissing_income,
    .groups = "drop"
  )

income_outlier_summary %>%
  mutate(
    income_p99_by_year = round(income_p99_by_year, 0),
    pct_removed = round(pct_removed, 2)
  ) %>%
  kable(
    caption = "Income outlier removal: observations above year-specific 99th percentile"
  )
Income outlier removal: observations above year-specific 99th percentile
year total_records records_with_nonmissing_income income_p99_by_year records_removed_as_income_outliers pct_removed
2019 1460 1240 33036 11 0.89
2021 2493 2293 25000 19 0.83
2023 313 313 64666 4 1.28
2024 711 711 394450 8 1.13
2025 1865 1865 54544 18 0.97
2026 455 455 33588 4 0.88

2. Real household income over time

Median real monthly household income by year

set.seed(123)

income_year <- hhs %>%
  filter(!is.na(real_monthly_hh_income_2021_mzn_clean)) %>%
  group_by(year) %>%
  summarise(
    n = n(),
    q25 = quantile(real_monthly_hh_income_2021_mzn_clean, 0.25, na.rm = TRUE),
    q75 = quantile(real_monthly_hh_income_2021_mzn_clean, 0.75, na.rm = TRUE),
    median_ci = list(bootstrap_median_ci(real_monthly_hh_income_2021_mzn_clean)),
    .groups = "drop"
  ) %>%
  unnest(median_ci)

income_year %>%
  mutate(
    across(c(q25, median, q75, ci_low, ci_high), ~ round(.x, 0))
  ) %>%
  kable(
    caption = "Median real monthly household income by year, constant 2021 MZN"
  )
Median real monthly household income by year, constant 2021 MZN
year n q25 q75 median ci_low ci_high
2019 1229 3304 9911 5506 5506 5506
2021 2274 1000 6500 3500 3500 3500
2023 309 7618 16930 11004 10158 12698
2024 703 3253 10573 4880 4880 4880
2025 1847 2338 7792 4675 4519 5065
2026 451 597 5971 1866 1493 2239 
ggplot(
  income_year,
  aes(x = year, y = median)
) +
  geom_ribbon(
    aes(ymin = ci_low, ymax = ci_high),
    alpha = 0.15
  ) +
  geom_line(linewidth = 0.9) +
  geom_point(aes(size = n), alpha = 0.85) +
  scale_x_continuous(breaks = sort(unique(income_year$year))) +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Median real monthly household income by year",
    subtitle = "Constant 2021 MZN; ribbon shows bootstrap 95% CI for the median; values above year-specific p99 excluded",
    x = "Survey year",
    y = "Median real monthly household income, 2021 MZN",
    size = "N"
  ) +
  theme(legend.position = "bottom")

Median real monthly household income by municipality and year

set.seed(123)

income_municipality_year <- hhs %>%
  filter(
    !is.na(g1_municipality),
    !is.na(year),
    !is.na(real_monthly_hh_income_2021_mzn_clean)
  ) %>%
  group_by(g1_municipality, year) %>%
  summarise(
    n = n(),
    median_ci = list(bootstrap_median_ci(real_monthly_hh_income_2021_mzn_clean)),
    .groups = "drop"
  ) %>%
  unnest(median_ci)

income_heatmap <- income_municipality_year %>%
  mutate(
    municipality = fct_reorder(g1_municipality, median, .fun = stats::median, na.rm = TRUE),
    label = paste0(comma(round(median, 0)), "\n", "n=", n),
    value_scaled = rescale(median, to = c(0, 1), from = range(median, na.rm = TRUE)),
    label_color = if_else(value_scaled < 0.45, "white", "black")
  )

ggplot(
  income_heatmap,
  aes(
    x = factor(year),
    y = municipality,
    fill = median
  )
) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 3, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(labels = comma, option = "C") +
  labs(
    title = "Median real monthly household income by municipality and year",
    subtitle = "Constant 2021 MZN; labels show median income and sample size",
    x = "Survey year",
    y = NULL,
    fill = "Median income"
  ) +
  theme(panel.grid = element_blank())

3. Income-source dependence and total household income

This section asks whether households with a higher share of income from artisanal fishing or farming tend to report higher or lower total real monthly household income. The analysis uses the cleaned real-income variable, so income is expressed in constant 2021 MZN and values above the year-specific 99th percentile are excluded.

The binned plots are the easiest to interpret descriptively. The adjusted models add controls for survey year and municipality, which helps reduce confounding from differences in where and when households were surveyed. However, these are still associations and should not be interpreted as causal effects of relying on fishing or farming.

Descriptive relationship using dependence bins

set.seed(123)

income_dependency_long <- hhs %>%
  filter(
    source_module_available,
    !is.na(real_monthly_hh_income_2021_mzn_clean)
  ) %>%
  transmute(
    year,
    g1_municipality,
    real_monthly_hh_income_2021_mzn_clean,
    `Artisanal fishing` = fishing_artisanal_income_pct,
    Farming = farming_income_pct
  ) %>%
  pivot_longer(
    cols = c(`Artisanal fishing`, Farming),
    names_to = "dependency_source",
    values_to = "income_dependency_pct"
  ) %>%
  mutate(
    dependency_bin = dependency_bins(income_dependency_pct)
  )

income_dependency_binned <- income_dependency_long %>%
  filter(!is.na(dependency_bin)) %>%
  group_by(dependency_source, dependency_bin) %>%
  summarise(
    n = n(),
    q25 = quantile(real_monthly_hh_income_2021_mzn_clean, 0.25, na.rm = TRUE),
    q75 = quantile(real_monthly_hh_income_2021_mzn_clean, 0.75, na.rm = TRUE),
    median_ci = list(bootstrap_median_ci(real_monthly_hh_income_2021_mzn_clean)),
    .groups = "drop"
  ) %>%
  unnest(median_ci)

income_dependency_binned %>%
  mutate(
    across(c(q25, median, q75, ci_low, ci_high), ~ round(.x, 0))
  ) %>%
  kable(
    caption = "Median real monthly household income by income-source dependence bin"
  )
Median real monthly household income by income-source dependence bin
dependency_source dependency_bin n q25 q75 median ci_low ci_high
Artisanal fishing 0% 2889 2000 7708 4500 4295 4675
Artisanal fishing 1–25% 471 3623 11008 7013 6772 7618
Artisanal fishing 26–50% 1133 1558 8133 4286 4066 4880
Artisanal fishing 51–75% 967 1558 7464 3800 3304 4000
Artisanal fishing 76–100% 1353 2500 10000 5065 5000 5506
Farming 0% 2980 3000 9653 5600 5506 6000
Farming 1–25% 1258 2903 10000 5000 4880 5506
Farming 26–50% 1686 1558 7000 3641 3500 3896
Farming 51–75% 466 653 5000 2120 1652 2338
Farming 76–100% 423 355 3747 1600 1493 2000 
ggplot(
  income_dependency_binned,
  aes(
    x = dependency_bin,
    y = median,
    group = dependency_source
  )
) +
  geom_errorbar(
    aes(ymin = ci_low, ymax = ci_high),
    width = 0.15,
    alpha = 0.55
  ) +
  geom_line(linewidth = 0.8) +
  geom_point(aes(size = n), alpha = 0.85) +
  facet_wrap(~ dependency_source, ncol = 1) +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Total real household income by income-source dependence",
    subtitle = "Points show median real monthly income; error bars show bootstrap 95% CIs for the median",
    x = "Share of household income from source",
    y = "Median real monthly household income, 2021 MZN",
    size = "N"
  ) +
  theme(legend.position = "bottom")

Scatterplots with smoothed relationship

The scatterplots below show the household-level relationship. The smooth line is useful for detecting non-linear patterns, but it should be interpreted carefully because repeated cross-sectional samples vary by year and municipality.

income_dependency_scatter <- income_dependency_long %>%
  filter(
    !is.na(income_dependency_pct),
    !is.na(real_monthly_hh_income_2021_mzn_clean)
  )

ggplot(
  income_dependency_scatter,
  aes(
    x = income_dependency_pct,
    y = real_monthly_hh_income_2021_mzn_clean
  )
) +
  geom_point(alpha = 0.12, size = 1) +
  geom_smooth(method = "loess", se = TRUE, linewidth = 0.9) +
  facet_wrap(~ dependency_source, ncol = 1) +
  scale_x_continuous(
    limits = c(0, 100),
    labels = label_percent(scale = 1)
  ) +
  scale_y_continuous(
    labels = comma,
    breaks = seq(0, 30000, 5000)
  ) +
  coord_cartesian(
    ylim = c(0, 30000)
  ) +
  labs(
    title = "Household income versus income-source dependence",
    subtitle = "Real monthly income in constant 2021 MZN; y-axis zoomed to 0–30,000 for readability",
    x = "Share of household income from source",
    y = "Real monthly household income, 2021 MZN"
  ) +
  theme_minimal(base_size = 12)

Adjusted association between income dependence and total income

The models below estimate the association between income-source dependence and total household income, adjusting for survey year and municipality. Because household income is skewed, the outcome is the log of real monthly household income. The coefficient is converted into an approximate percent difference in total income for each 10 percentage point increase in income dependence.

income_dependency_model_data <- hhs %>%
  filter(
    source_module_available,
    !is.na(real_monthly_hh_income_2021_mzn_clean),
    real_monthly_hh_income_2021_mzn_clean >= 0,
    !is.na(g1_municipality),
    !is.na(year)
  ) %>%
  mutate(
    log_real_income = log1p(real_monthly_hh_income_2021_mzn_clean),
    fishing_dependency_10pp = fishing_artisanal_income_pct / 10,
    farming_dependency_10pp = farming_income_pct / 10,
    year_factor = factor(year),
    municipality_factor = factor(g1_municipality)
  )

fit_income_dependency_model <- function(data, dependency_var, source_label) {

  model_formula <- as.formula(
    paste0("log_real_income ~ ", dependency_var, " + year_factor + municipality_factor")
  )

  model <- lm(model_formula, data = data)

  tidy(model, conf.int = TRUE) %>%
    filter(term == dependency_var) %>%
    transmute(
      dependency_source = source_label,
      percent_difference_per_10pp = 100 * (exp(estimate) - 1),
      ci_low = 100 * (exp(conf.low) - 1),
      ci_high = 100 * (exp(conf.high) - 1),
      p_value = p.value
    )
}

income_dependency_model_results <- bind_rows(
  fit_income_dependency_model(
    income_dependency_model_data,
    dependency_var = "fishing_dependency_10pp",
    source_label = "Artisanal fishing"
  ),
  fit_income_dependency_model(
    income_dependency_model_data,
    dependency_var = "farming_dependency_10pp",
    source_label = "Farming"
  )
)

income_dependency_model_results %>%
  mutate(
    percent_difference_per_10pp = round(percent_difference_per_10pp, 1),
    ci_low = round(ci_low, 1),
    ci_high = round(ci_high, 1),
    p_value = scales::pvalue(p_value)
  ) %>%
  kable(
    caption = "Adjusted association between income-source dependence and total real household income"
  )
Adjusted association between income-source dependence and total real household income
dependency_source percent_difference_per_10pp ci_low ci_high p_value
Artisanal fishing 3.1 2.1 4 <0.001
Farming -11.2 -12.4 -10 <0.001 
ggplot(
  income_dependency_model_results,
  aes(
    x = percent_difference_per_10pp,
    y = dependency_source,
    xmin = ci_low,
    xmax = ci_high
  )
) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.6) +
  geom_errorbarh(height = 0.15, alpha = 0.7) +
  geom_point(size = 2.8) +
  labs(
    title = "Adjusted association: income dependence and total household income",
    subtitle = "Percent difference in real monthly household income per 10 percentage point increase; adjusted for year and municipality",
    x = "% difference in total real monthly household income",
    y = NULL
  )

Interpretation note

Households that rely more on artisanal fishing tend to have slightly higher total real monthly household income, while households that rely more on farming tend to have much lower total real monthly household income, after adjusting for year and municipality.

More specifically:

  • For artisanal fishing, every 10 percentage-point increase in the share of household income coming from artisanal fishing is associated with about 3.1% higher total real monthly household income. The confidence interval is roughly 2.1% to 4.0%, and the p-value is <0.001, so this is statistically significant.

  • For farming, every 10 percentage-point increase in the share of household income coming from farming is associated with about 11.2% lower total real monthly household income. The confidence interval is roughly -12.4% to -10.0%, and the p-value is <0.001, so also statistically significant.

4. Income-source composition over time

Whole Mozambique: mean reported income share by source and year

source_long <- hhs %>%
  filter(source_module_available) %>%
  select(
    year,
    g1_province,
    g1_municipality,
    g1_community,
    all_of(source_cols),
    food_secure,
    income_sufficient
  ) %>%
  pivot_longer(
    cols = all_of(source_cols),
    names_to = "source_col",
    values_to = "income_share_pct"
  ) %>%
  left_join(source_labels, by = "source_col")

source_year <- source_long %>%
  group_by(year, source) %>%
  summarise(
    n = n(),
    mean_share = mean(income_share_pct, na.rm = TRUE),
    median_share = median(income_share_pct, na.rm = TRUE),
    pct_households_any_income = 100 * mean(income_share_pct > 0, na.rm = TRUE),
    .groups = "drop"
  )

source_year %>%
  mutate(
    mean_share = round(mean_share, 1),
    median_share = round(median_share, 1),
    pct_households_any_income = round(pct_households_any_income, 1)
  ) %>%
  arrange(year, desc(mean_share)) %>%
  kable(
    caption = "Income-source summary by year"
  )
Income-source summary by year
year source n mean_share median_share pct_households_any_income
2019 Artisanal fishing 1460 32.7 0 41.3
2019 Farming 1460 18.2 0 40.1
2019 Other 1460 15.1 0 22.1
2019 Other wage 1460 6.3 0 8.9
2019 Aquaculture 1460 2.6 0 5.9
2019 Tourism 1460 2.0 0 2.9
2019 Harvesting 1460 0.6 0 1.4
2019 Extraction 1460 0.4 0 1.5
2019 Buying/trading 1460 0.0 0 0.0
2019 Industrial 1460 0.0 0 0.0
2019 Processing 1460 0.0 0 0.0
2021 Artisanal fishing 2493 34.8 0 45.8
2021 Buying/trading 2493 29.0 0 38.9
2021 Farming 2493 22.2 10 52.1
2021 Other 2493 3.5 0 8.7
2021 Processing 2493 2.2 0 5.0
2021 Aquaculture 2493 1.6 0 5.7
2021 Harvesting 2493 1.2 0 3.9
2021 Other wage 2493 0.3 0 0.5
2021 Industrial 2493 0.2 0 0.9
2021 Tourism 2493 0.2 0 1.0
2021 Extraction 2493 0.0 0 0.0
2023 Artisanal fishing 313 33.5 25 82.7
2023 Farming 313 21.7 25 60.4
2023 Industrial 313 14.3 0 41.5
2023 Harvesting 313 10.0 0 30.4
2023 Buying/trading 313 9.8 0 35.5
2023 Processing 313 6.0 0 26.2
2023 Tourism 313 1.8 0 9.9
2023 Aquaculture 313 1.4 0 6.1
2023 Other wage 313 0.9 0 8.0
2023 Extraction 313 0.5 0 2.9
2023 Other 313 0.1 0 0.6
2024 Artisanal fishing 711 44.7 50 78.2
2024 Farming 711 28.4 25 84.0
2024 Buying/trading 711 18.6 15 65.5
2024 Processing 711 7.1 0 31.9
2024 Aquaculture 711 0.6 0 1.8
2024 Harvesting 711 0.3 0 1.8
2024 Other 711 0.2 0 0.8
2024 Extraction 711 0.0 0 0.0
2024 Industrial 711 0.0 0 0.1
2024 Other wage 711 0.0 0 0.1
2024 Tourism 711 0.0 0 0.3
2025 Artisanal fishing 1865 39.1 35 67.0
2025 Farming 1865 20.7 15 52.4
2025 Buying/trading 1865 20.2 0 43.9
2025 Processing 1865 7.6 0 23.7
2025 Other 1865 4.9 0 9.0
2025 Industrial 1865 2.2 0 9.8
2025 Aquaculture 1865 1.8 0 7.8
2025 Other wage 1865 1.6 0 2.7
2025 Harvesting 1865 1.2 0 5.1
2025 Extraction 1865 0.4 0 1.6
2025 Tourism 1865 0.3 0 0.8
2026 Farming 455 43.4 40 85.1
2026 Artisanal fishing 455 39.2 40 74.9
2026 Buying/trading 455 10.3 0 16.3
2026 Other 455 4.1 0 11.4
2026 Harvesting 455 2.2 0 7.0
2026 Processing 455 0.4 0 1.1
2026 Other wage 455 0.3 0 0.7
2026 Aquaculture 455 0.2 0 0.7
2026 Extraction 455 0.0 0 0.0
2026 Industrial 455 0.0 0 0.2
2026 Tourism 455 0.0 0 0.0 
source_year_heatmap <- source_year %>%
  mutate(
    source = fct_reorder(source, mean_share, .fun = mean, na.rm = TRUE),
    label = round(mean_share, 1),
    value_scaled = rescale(mean_share, to = c(0, 1), from = range(mean_share, na.rm = TRUE)),
    label_color = if_else(value_scaled < 0.45, "white", "black")
  )

ggplot(
  source_year_heatmap,
  aes(
    x = factor(year),
    y = source,
    fill = mean_share
  )
) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 3, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(option = "C", labels = label_percent(scale = 1)) +
  labs(
    title = "Mean reported share of household income by source and year",
    subtitle = "Values are percentages of total household income reported for each source",
    x = "Survey year",
    y = NULL,
    fill = "Mean share"
  ) +
  theme(panel.grid = element_blank())

top_sources <- source_long %>%
  group_by(source) %>%
  summarise(overall_mean_share = mean(income_share_pct, na.rm = TRUE), .groups = "drop") %>%
  slice_max(overall_mean_share, n = 6) %>%
  pull(source)

source_year %>%
  filter(source %in% top_sources) %>%
  ggplot(
    aes(
      x = year,
      y = mean_share,
      color = source,
      group = source
    )
  ) +
  geom_line(linewidth = 0.9) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = sort(unique(source_year$year))) +
  scale_y_continuous(labels = label_percent(scale = 1)) +
  scale_color_brewer(palette = "Dark2") +
  labs(
    title = "Top income-source trends across Mozambique",
    subtitle = "Mean reported share of household income; top six sources by overall mean share",
    x = "Survey year",
    y = "Mean share of household income",
    color = "Income source"
  ) +
  theme(legend.position = "bottom")

By municipality and year: selected income sources

The following heatmaps focus on the main income sources most relevant to the analysis. They show the mean reported share of income from each source by municipality-year.

source_municipality_year <- source_long %>%
  group_by(g1_municipality, year, source) %>%
  summarise(
    n = n(),
    mean_share = mean(income_share_pct, na.rm = TRUE),
    pct_households_any_income = 100 * mean(income_share_pct > 0, na.rm = TRUE),
    .groups = "drop"
  )

plot_source_municipality_year <- function(source_name) {

  plot_data <- source_municipality_year %>%
    filter(source == source_name) %>%
    mutate(
      municipality = fct_reorder(g1_municipality, mean_share, .fun = mean, na.rm = TRUE),
      label = paste0(round(mean_share, 1), "\n", "n=", n),
      value_scaled = rescale(mean_share, to = c(0, 1), from = range(mean_share, na.rm = TRUE)),
      label_color = if_else(value_scaled < 0.45, "white", "black")
    )

  ggplot(
    plot_data,
    aes(
      x = factor(year),
      y = municipality,
      fill = mean_share
    )
  ) +
    geom_tile(color = "white") +
    geom_text(aes(label = label, color = label_color), size = 3, fontface = "bold") +
    scale_color_identity() +
    scale_fill_viridis_c(option = "C", labels = label_percent(scale = 1)) +
    labs(
      title = paste0("Mean income share from ", source_name, " by municipality and year"),
      subtitle = "Cell labels show mean income share (%) and HHS records",
      x = "Survey year",
      y = NULL,
      fill = "Mean share"
    ) +
    theme(panel.grid = element_blank())
}
plot_source_municipality_year("Artisanal fishing")

plot_source_municipality_year("Farming")

plot_source_municipality_year("Buying/trading")

5. Artisanal fishing income dependence

Whole Mozambique by year

fishing_year <- hhs %>%
  filter(source_module_available) %>%
  group_by(year) %>%
  summarise(
    n = n(),
    mean_ci = list(mean_ci(fishing_artisanal_income_pct)),
    median_ci = list(bootstrap_median_ci(fishing_artisanal_income_pct)),
    pct_any_fishing = 100 * mean(fishing_artisanal_income_pct > 0, na.rm = TRUE),
    pct_high_fishing_50 = 100 * mean(fishing_artisanal_income_pct >= 50, na.rm = TRUE),
    pct_very_high_fishing_75 = 100 * mean(fishing_artisanal_income_pct >= 75, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  unnest_wider(mean_ci, names_sep = "_") %>%
  unnest_wider(median_ci, names_sep = "_") %>%
  rename(
    mean_fishing_share = mean_ci_mean,
    mean_ci_low = mean_ci_ci_low,
    mean_ci_high = mean_ci_ci_high,
    median_fishing_share = median_ci_median,
    median_ci_low = median_ci_ci_low,
    median_ci_high = median_ci_ci_high
  )

fishing_year %>%
  mutate(
    across(
      c(
        mean_fishing_share,
        mean_ci_low,
        mean_ci_high,
        median_fishing_share,
        median_ci_low,
        median_ci_high,
        pct_any_fishing,
        pct_high_fishing_50,
        pct_very_high_fishing_75
      ),
      ~ round(.x, 1)
    )
  ) %>%
  kable(
    caption = "Artisanal fishing income dependence by year"
  )
Artisanal fishing income dependence by year
year n mean_ci_n mean_fishing_share mean_ci_low mean_ci_high median_fishing_share median_ci_low median_ci_high pct_any_fishing pct_high_fishing_50 pct_very_high_fishing_75
2019 1460 1460 32.7 30.5 34.8 0 0 0 41.3 37.3 27.8
2021 2493 2493 34.8 33.1 36.4 0 0 0 45.8 38.5 26.8
2023 313 313 33.5 31.1 35.9 25 25 30 82.7 41.9 8.9
2024 711 711 44.7 42.5 46.9 50 50 50 78.2 61.6 16.2
2025 1865 1865 39.1 37.6 40.7 35 30 40 67.0 44.2 19.1
2026 455 455 39.2 36.4 41.9 40 40 40 74.9 40.7 12.7 
ggplot(
  fishing_year,
  aes(x = year, y = mean_fishing_share)
) +
  geom_ribbon(aes(ymin = mean_ci_low, ymax = mean_ci_high), alpha = 0.15) +
  geom_line(linewidth = 0.9) +
  geom_point(aes(size = n), alpha = 0.85) +
  scale_x_continuous(breaks = sort(unique(fishing_year$year))) +
  scale_y_continuous(limits = c(0, 100), labels = label_percent(scale = 1)) +
  labs(
    title = "Mean share of household income from artisanal fishing by year",
    subtitle = "Ribbon shows approximate 95% CI for the mean",
    x = "Survey year",
    y = "Mean share of income from artisanal fishing",
    size = "N"
  ) +
  theme(legend.position = "bottom")

Artisanal fishing dependence categories

fishing_dependence_categories <- hhs %>%
  filter(source_module_available) %>%
  mutate(
    fishing_dependence_bin = dependency_bins(fishing_artisanal_income_pct)
  ) %>%
  count(year, fishing_dependence_bin, name = "n") %>%
  group_by(year) %>%
  mutate(pct = 100 * n / sum(n)) %>%
  ungroup()

ggplot(
  fishing_dependence_categories,
  aes(
    x = factor(year),
    y = pct,
    fill = fishing_dependence_bin
  )
) +
  geom_col(position = "fill") +
  scale_y_continuous(labels = percent) +
  scale_fill_brewer(palette = "Blues", direction = 1) +
  labs(
    title = "Distribution of artisanal fishing income dependence by year",
    x = "Survey year",
    y = "Share of households",
    fill = "Income from artisanal fishing"
  ) +
  theme(legend.position = "bottom")

Artisanal fishing dependence by municipality and year

fishing_municipality_year <- hhs %>%
  filter(
    source_module_available,
    !is.na(g1_municipality),
    !is.na(year)
  ) %>%
  group_by(g1_municipality, year) %>%
  summarise(
    n = n(),
    mean_fishing_share = mean(fishing_artisanal_income_pct, na.rm = TRUE),
    pct_any_fishing = 100 * mean(fishing_artisanal_income_pct > 0, na.rm = TRUE),
    pct_high_fishing_50 = 100 * mean(fishing_artisanal_income_pct >= 50, na.rm = TRUE),
    .groups = "drop"
  )

fishing_mun_heatmap <- fishing_municipality_year %>%
  mutate(
    municipality = fct_reorder(g1_municipality, mean_fishing_share, .fun = mean, na.rm = TRUE),
    label = paste0(round(mean_fishing_share, 1), "\n", "n=", n),
    value_scaled = rescale(mean_fishing_share, to = c(0, 1), from = range(mean_fishing_share, na.rm = TRUE)),
    label_color = if_else(value_scaled < 0.45, "white", "black")
  )

ggplot(
  fishing_mun_heatmap,
  aes(
    x = factor(year),
    y = municipality,
    fill = mean_fishing_share
  )
) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 3, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(option = "C", labels = label_percent(scale = 1)) +
  labs(
    title = "Mean share of income from artisanal fishing by municipality and year",
    subtitle = "Cell labels show mean income share (%) and HHS records",
    x = "Survey year",
    y = NULL,
    fill = "Mean share"
  ) +
  theme(panel.grid = element_blank())

6. Farming income dependence

This section mirrors the artisanal fishing analysis, but for the share of income reported from farming.

farming_year <- hhs %>%
  filter(source_module_available) %>%
  group_by(year) %>%
  summarise(
    n = n(),
    mean_ci = list(mean_ci(farming_income_pct)),
    median_ci = list(bootstrap_median_ci(farming_income_pct)),
    pct_any_farming = 100 * mean(farming_income_pct > 0, na.rm = TRUE),
    pct_high_farming_50 = 100 * mean(farming_income_pct >= 50, na.rm = TRUE),
    pct_very_high_farming_75 = 100 * mean(farming_income_pct >= 75, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  unnest_wider(mean_ci, names_sep = "_") %>%
  unnest_wider(median_ci, names_sep = "_") %>%
  rename(
    mean_farming_share = mean_ci_mean,
    mean_ci_low = mean_ci_ci_low,
    mean_ci_high = mean_ci_ci_high,
    median_farming_share = median_ci_median,
    median_ci_low = median_ci_ci_low,
    median_ci_high = median_ci_ci_high
  )

farming_year %>%
  mutate(
    across(
      c(
        mean_farming_share,
        mean_ci_low,
        mean_ci_high,
        median_farming_share,
        median_ci_low,
        median_ci_high,
        pct_any_farming,
        pct_high_farming_50,
        pct_very_high_farming_75
      ),
      ~ round(.x, 1)
    )
  ) %>%
  kable(
    caption = "Farming income dependence by year"
  )
Farming income dependence by year
year n mean_ci_n mean_farming_share mean_ci_low mean_ci_high median_farming_share median_ci_low median_ci_high pct_any_farming pct_high_farming_50 pct_very_high_farming_75
2019 1460 1460 18.2 16.8 19.7 0 0.0 0 40.1 15.5 7.9
2021 2493 2493 22.2 21.1 23.4 10 1.0 10 52.1 18.4 8.5
2023 313 313 21.7 19.4 24.1 25 24.9 25 60.4 22.7 1.0
2024 711 711 28.4 26.9 29.9 25 25.0 25 84.0 17.9 4.2
2025 1865 1865 20.7 19.5 21.8 15 5.0 20 52.4 16.1 4.0
2026 455 455 43.4 40.8 46.0 40 40.0 40 85.1 44.4 13.0 
ggplot(
  farming_year,
  aes(x = year, y = mean_farming_share)
) +
  geom_ribbon(aes(ymin = mean_ci_low, ymax = mean_ci_high), alpha = 0.15) +
  geom_line(linewidth = 0.9) +
  geom_point(aes(size = n), alpha = 0.85) +
  scale_x_continuous(breaks = sort(unique(farming_year$year))) +
  scale_y_continuous(limits = c(0, 100), labels = label_percent(scale = 1)) +
  labs(
    title = "Mean share of household income from farming by year",
    subtitle = "Ribbon shows approximate 95% CI for the mean",
    x = "Survey year",
    y = "Mean share of income from farming",
    size = "N"
  ) +
  theme(legend.position = "bottom")

farming_municipality_year <- hhs %>%
  filter(
    source_module_available,
    !is.na(g1_municipality),
    !is.na(year)
  ) %>%
  group_by(g1_municipality, year) %>%
  summarise(
    n = n(),
    mean_farming_share = mean(farming_income_pct, na.rm = TRUE),
    pct_any_farming = 100 * mean(farming_income_pct > 0, na.rm = TRUE),
    pct_high_farming_50 = 100 * mean(farming_income_pct >= 50, na.rm = TRUE),
    .groups = "drop"
  )

farming_mun_heatmap <- farming_municipality_year %>%
  mutate(
    municipality = fct_reorder(g1_municipality, mean_farming_share, .fun = mean, na.rm = TRUE),
    label = paste0(round(mean_farming_share, 1), "\n", "n=", n),
    value_scaled = rescale(mean_farming_share, to = c(0, 1), from = range(mean_farming_share, na.rm = TRUE)),
    label_color = if_else(value_scaled < 0.45, "white", "black")
  )

ggplot(
  farming_mun_heatmap,
  aes(
    x = factor(year),
    y = municipality,
    fill = mean_farming_share
  )
) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 3, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(option = "C", labels = label_percent(scale = 1)) +
  labs(
    title = "Mean share of income from farming by municipality and year",
    subtitle = "Cell labels show mean income share (%) and HHS records",
    x = "Survey year",
    y = NULL,
    fill = "Mean share"
  ) +
  theme(panel.grid = element_blank())

7. Relationship between income dependence and Sustainable Livelihoods

Descriptive binned relationships

The figures below show how the two Sustainable Livelihoods indicators vary across income-dependence bins. The bins are descriptive only and do not adjust for differences in year, municipality, or sample composition.

relationship_binned <- hhs %>%
  filter(source_module_available) %>%
  select(
    year,
    g1_municipality,
    fishing_artisanal_income_pct,
    farming_income_pct,
    food_secure,
    income_sufficient
  ) %>%
  mutate(
    `Artisanal fishing` = fishing_artisanal_income_pct,
    Farming = farming_income_pct
  ) %>%
  pivot_longer(
    cols = c(`Artisanal fishing`, Farming),
    names_to = "dependency_source",
    values_to = "income_dependency_pct"
  ) %>%
  mutate(
    dependency_bin = dependency_bins(income_dependency_pct)
  ) %>%
  pivot_longer(
    cols = c(food_secure, income_sufficient),
    names_to = "outcome",
    values_to = "positive"
  ) %>%
  mutate(
    outcome = recode(
      outcome,
      food_secure = "Food security",
      income_sufficient = "Income sufficiency / ability to meet needs"
    )
  ) %>%
  group_by(dependency_source, dependency_bin, outcome) %>%
  summarise(
    prop = list(prop_ci(positive)),
    .groups = "drop"
  ) %>%
  unnest(prop)

relationship_binned %>%
  mutate(
    pct = round(pct, 1),
    ci_low = round(ci_low, 1),
    ci_high = round(ci_high, 1)
  ) %>%
  kable(
    caption = "Sustainable Livelihoods indicators by income-dependence bin"
  )
Sustainable Livelihoods indicators by income-dependence bin
dependency_source dependency_bin outcome n pct ci_low ci_high
Artisanal fishing 0% Food security 3081 12.3 11.2 13.5
Artisanal fishing 0% Income sufficiency / ability to meet needs 3079 28.6 27.0 30.2
Artisanal fishing 1–25% Food security 499 7.2 4.9 9.5
Artisanal fishing 1–25% Income sufficiency / ability to meet needs 498 25.1 21.3 28.9
Artisanal fishing 26–50% Food security 1159 10.2 8.4 11.9
Artisanal fishing 26–50% Income sufficiency / ability to meet needs 1156 19.1 16.9 21.4
Artisanal fishing 51–75% Food security 997 10.8 8.9 12.8
Artisanal fishing 51–75% Income sufficiency / ability to meet needs 995 22.3 19.7 24.9
Artisanal fishing 76–100% Food security 1470 13.9 12.1 15.6
Artisanal fishing 76–100% Income sufficiency / ability to meet needs 1445 31.4 29.0 33.8
Farming 0% Food security 3215 10.7 9.6 11.8
Farming 0% Income sufficiency / ability to meet needs 3178 28.9 27.3 30.4
Farming 1–25% Food security 1312 16.6 14.6 18.6
Farming 1–25% Income sufficiency / ability to meet needs 1312 34.8 32.2 37.3
Farming 26–50% Food security 1765 11.2 9.7 12.7
Farming 26–50% Income sufficiency / ability to meet needs 1756 18.5 16.6 20.3
Farming 51–75% Food security 480 10.0 7.3 12.7
Farming 51–75% Income sufficiency / ability to meet needs 481 28.3 24.2 32.3
Farming 76–100% Food security 434 8.8 6.1 11.4
Farming 76–100% Income sufficiency / ability to meet needs 446 15.7 12.3 19.1 
ggplot(
  relationship_binned,
  aes(
    x = dependency_bin,
    y = pct,
    color = outcome,
    group = outcome
  )
) +
  geom_errorbar(
    aes(ymin = ci_low, ymax = ci_high),
    width = 0.15,
    position = position_dodge(width = 0.25),
    alpha = 0.5
  ) +
  geom_line(position = position_dodge(width = 0.25), linewidth = 0.8) +
  geom_point(position = position_dodge(width = 0.25), size = 2) +
  facet_wrap(~ dependency_source, ncol = 1) +
  scale_y_continuous(limits = c(0, 100), labels = label_percent(scale = 1)) +
  scale_color_brewer(palette = "Dark2") +
  labs(
    title = "Sustainable Livelihoods indicators by income-source dependence",
    subtitle = "Error bars show approximate 95% CIs for the proportion",
    x = "Share of household income from source",
    y = "% positive",
    color = "Indicator"
  ) +
  theme(legend.position = "bottom")

Exploratory adjusted associations

The models below estimate the association between income-source dependence and each Sustainable Livelihoods indicator, adjusting for survey year and municipality. The coefficient is shown as an odds ratio for each 10 percentage point increase in income dependence.

These models are exploratory. They should not be interpreted causally because income dependence, livelihood outcomes, and site/year conditions are likely jointly determined.

model_data <- hhs %>%
  filter(source_module_available) %>%
  mutate(
    fishing_dependency_10pp = fishing_artisanal_income_pct / 10,
    farming_dependency_10pp = farming_income_pct / 10,
    year_factor = factor(year),
    municipality_factor = factor(g1_municipality)
  )

fit_dependency_model <- function(data, outcome_var, dependency_var, source_label, outcome_label) {

  model_formula <- as.formula(
    paste0(outcome_var, " ~ ", dependency_var, " + year_factor + municipality_factor")
  )

  model <- glm(
    model_formula,
    data = data,
    family = binomial()
  )

  tidy(model, conf.int = TRUE, exponentiate = TRUE) %>%
    filter(term == dependency_var) %>%
    transmute(
      dependency_source = source_label,
      outcome = outcome_label,
      odds_ratio_per_10pp = estimate,
      ci_low = conf.low,
      ci_high = conf.high,
      p_value = p.value
    )
}

model_results <- bind_rows(
  fit_dependency_model(
    model_data,
    outcome_var = "food_secure",
    dependency_var = "fishing_dependency_10pp",
    source_label = "Artisanal fishing",
    outcome_label = "Food security"
  ),
  fit_dependency_model(
    model_data,
    outcome_var = "income_sufficient",
    dependency_var = "fishing_dependency_10pp",
    source_label = "Artisanal fishing",
    outcome_label = "Income sufficiency / ability to meet needs"
  ),
  fit_dependency_model(
    model_data,
    outcome_var = "food_secure",
    dependency_var = "farming_dependency_10pp",
    source_label = "Farming",
    outcome_label = "Food security"
  ),
  fit_dependency_model(
    model_data,
    outcome_var = "income_sufficient",
    dependency_var = "farming_dependency_10pp",
    source_label = "Farming",
    outcome_label = "Income sufficiency / ability to meet needs"
  )
)

model_results %>%
  mutate(
    odds_ratio_per_10pp = round(odds_ratio_per_10pp, 3),
    ci_low = round(ci_low, 3),
    ci_high = round(ci_high, 3),
    p_value = scales::pvalue(p_value)
  ) %>%
  kable(
    caption = "Exploratory adjusted associations: odds ratio per 10 percentage point increase in income dependence"
  )
Exploratory adjusted associations: odds ratio per 10 percentage point increase in income dependence
dependency_source outcome odds_ratio_per_10pp ci_low ci_high p_value
Artisanal fishing Food security 0.985 0.964 1.005 0.140
Artisanal fishing Income sufficiency / ability to meet needs 0.988 0.973 1.002 0.101
Farming Food security 0.918 0.888 0.948 <0.001
Farming Income sufficiency / ability to meet needs 0.874 0.853 0.895 <0.001 
ggplot(
  model_results,
  aes(
    x = odds_ratio_per_10pp,
    y = outcome,
    xmin = ci_low,
    xmax = ci_high
  )
) +
  geom_vline(xintercept = 1, linetype = "dashed", alpha = 0.6) +
  geom_errorbarh(height = 0.15, alpha = 0.7) +
  geom_point(size = 2.5) +
  facet_wrap(~ dependency_source, ncol = 1) +
  scale_x_log10() +
  labs(
    title = "Adjusted association between income-source dependence and Sustainable Livelihoods",
    subtitle = "Odds ratios per 10 percentage point increase; adjusted for year and municipality",
    x = "Odds ratio, log scale",
    y = NULL
  )

Interpretation note

Greater farming income dependence is associated with worse Sustainable Livelihoods outcomes, while artisanal fishing dependence does not show a statistically clear relationship with those outcomes after adjusting for year and municipality.

More specifically:

  • For artisanal fishing, the odds ratios are very close to 1:

Food security: OR = 0.985, p = 0.140 Income sufficiency / ability to meet needs: OR = 0.988, p = 0.101

This means that a 10 percentage-point increase in the share of income from artisanal fishing is associated with slightly lower odds of positive livelihood outcomes, but the confidence intervals cross 1 and the p-values are not statistically significant (i.e., no clear evidence of an association).

  • For farming, the relationship is stronger and statistically significant:

Food security: OR = 0.918, p < 0.001 Income sufficiency / ability to meet needs: OR = 0.874, p < 0.001

This means that for every 10 percentage-point increase in the share of income coming from farming, households have about 8.2% lower odds of reporting food security and about 12.6% lower odds of reporting that they can meet household needs, after adjusting for year and municipality.

NOTE: The descriptive binned plot does not show a strong monotonic decline because it pools observations across municipalities and years and summarizes a continuous variable into broad categories. However, the adjusted model compares households after accounting for year and municipality, and in that specification greater farming dependence is consistently associated with lower odds of food security and income sufficiency. This clearly suggests that the farming-dependence relationship is strong within comparable municipality-year contexts.

8. Sustainable Livelihoods through time: are negative conditions becoming less severe?

The CCRF Sustainable Livelihoods score is intentionally strict: households are counted as positive only if they report never worrying about food and if they can meet household needs fairly easily, easily, or very easily.

However, for adaptive management it is also useful to ask a more diagnostic question: are the underlying conditions becoming less negative through time? For this section, the analysis therefore looks directly at two negative indicators:

  1. Food worry: the household reports that worrying about not having enough food was Sometimes or Often true during the last 12 months.
  2. Financial strain: the household reports covering needs With difficulty or With great difficulty.

A decline in these negative indicators over time would suggest improvement, even if the strict positive Sustainable Livelihoods score remains low. The results should still be interpreted cautiously because the HHS is repeated cross-sectional data and the set of municipalities/sites sampled changes across years.

First-to-latest observed change by municipality

This plot summarizes whether each municipality looks more or less negative between the first and latest HHS year available for that municipality. Because municipalities are not observed in every year, this is a simple descriptive comparison, not a formal trend estimate.

sl_negative_change <- sl_municipality_year %>%
  filter(indicator_direction == "Negative indicator") %>%
  filter(!is.na(pct)) %>%
  arrange(g1_municipality, indicator, year) %>%
  group_by(g1_municipality, indicator) %>%
  summarise(
    first_year = first(year),
    last_year = last(year),
    first_pct = first(pct),
    last_pct = last(pct),
    first_n = first(n),
    last_n = last(n),
    n_years = n_distinct(year),
    change_pp = last_pct - first_pct,
    .groups = "drop"
  ) %>%
  filter(n_years >= 2) %>%
  mutate(
    change_direction = case_when(
      change_pp < 0 ~ "Less negative",
      change_pp > 0 ~ "More negative",
      TRUE ~ "No change"
    ),
    change_label = paste0(round(change_pp, 1), " pp")
  )

sl_negative_change %>%
  arrange(indicator, change_pp) %>%
  mutate(
    first_pct = round(first_pct, 1),
    last_pct = round(last_pct, 1),
    change_pp = round(change_pp, 1)
  ) %>%
  kable(caption = "Change in negative Sustainable Livelihoods indicators from first to latest observed year")
Change in negative Sustainable Livelihoods indicators from first to latest observed year
g1_municipality indicator first_year last_year first_pct last_pct first_n last_n n_years change_pp change_direction change_label
Inharrime Food worry 2019 2021 78.4 77.8 185 225 2 -0.6 Less negative -0.6 pp
Dondo Food worry 2019 2021 100.0 100.0 199 27 2 0.0 No change 0 pp
Inhassoro Food worry 2019 2025 90.2 90.7 194 600 3 0.5 More negative 0.5 pp
Memba Food worry 2019 2026 95.0 97.4 201 274 4 2.4 More negative 2.4 pp
Nacala Porto Food worry 2024 2026 80.0 86.2 260 181 2 6.2 More negative 6.2 pp
Ilha de Mocambique Food worry 2019 2025 89.6 99.3 327 957 4 9.7 More negative 9.7 pp
Massinga Food worry 2019 2021 70.6 93.4 136 137 2 22.8 More negative 22.8 pp
Mogincual Food worry 2024 2025 63.1 99.4 306 308 2 36.3 More negative 36.3 pp
Matutuíne Food worry 2019 2021 25.3 74.2 158 151 2 48.9 More negative 48.9 pp
Nacala Porto Financial strain 2024 2026 81.5 43.1 260 181 2 -38.4 Less negative -38.4 pp
Dondo Financial strain 2019 2021 90.4 53.3 198 15 2 -37.1 Less negative -37.1 pp
Mogincual Financial strain 2024 2025 69.6 56.2 306 308 2 -13.4 Less negative -13.4 pp
Massinga Financial strain 2019 2021 63.2 54.7 133 137 2 -8.4 Less negative -8.4 pp
Matutuíne Financial strain 2019 2021 56.6 53.3 152 150 2 -3.2 Less negative -3.2 pp
Inharrime Financial strain 2019 2021 58.3 56.6 204 219 2 -1.7 Less negative -1.7 pp
Inhassoro Financial strain 2019 2025 70.5 76.7 190 600 3 6.1 More negative 6.1 pp
Memba Financial strain 2019 2026 66.8 76.6 202 274 4 9.8 More negative 9.8 pp
Ilha de Mocambique Financial strain 2019 2025 59.1 83.2 325 957 4 24.1 More negative 24.1 pp 
ggplot(
  sl_negative_change,
  aes(
    x = change_pp,
    y = fct_reorder(g1_municipality, change_pp),
    fill = change_direction
  )
) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey40") +
  geom_col(width = 0.7) +
  geom_text(
    aes(
      label = change_label,
      hjust = if_else(change_pp < 0, 1.1, -0.1)
    ),
    size = 3,
    show.legend = FALSE
  ) +
  facet_wrap(~ indicator, ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = function(x) paste0(x, " pp")) +
  labs(
    title = "Change in negative Sustainable Livelihoods conditions",
    subtitle = "Negative values mean the condition became less common between the first and latest observed HHS year",
    x = "Change in percentage points",
    y = NULL,
    fill = "Direction"
  ) +
  theme(legend.position = "bottom")

Exploratory municipality-specific time-trend models

As a simple statistical check, the models below fit a separate logistic regression for each municipality and negative outcome:

\[ \Pr(\text{negative condition}=1) = f(\text{survey year}) \]

The reported odds ratio is the multiplicative change in the odds of the negative condition for each additional calendar year. Values below 1 suggest that the condition is becoming less negative through time; values above 1 suggest it is becoming more negative.

These models are exploratory and can be unstable where a municipality has few survey years, small sample sizes, or little variation in the outcome.

fit_municipality_trend <- function(data, outcome_var, outcome_label) {

  data %>%
    filter(
      !is.na(.data[[outcome_var]]),
      !is.na(year),
      !is.na(g1_municipality)
    ) %>%
    group_by(g1_municipality) %>%
    group_modify(function(.x, .y) {

      dd <- .x %>% mutate(year_numeric = as.numeric(year))

      if (
        nrow(dd) < 30 ||
        n_distinct(dd$year_numeric) < 2 ||
        n_distinct(dd[[outcome_var]]) < 2
      ) {
        return(
          tibble(
            outcome = outcome_label,
            n = nrow(dd),
            n_years = n_distinct(dd$year_numeric),
            odds_ratio_per_year = NA_real_,
            ci_low = NA_real_,
            ci_high = NA_real_,
            p_value = NA_real_,
            trend_interpretation = "Insufficient variation/data"
          )
        )
      }

      fit <- tryCatch(
        suppressWarnings(glm(as.formula(paste0(outcome_var, " ~ year_numeric")), data = dd, family = binomial())),
        error = function(e) NULL
      )

      if (is.null(fit) || !("year_numeric" %in% rownames(summary(fit)$coefficients))) {
        return(
          tibble(
            outcome = outcome_label,
            n = nrow(dd),
            n_years = n_distinct(dd$year_numeric),
            odds_ratio_per_year = NA_real_,
            ci_low = NA_real_,
            ci_high = NA_real_,
            p_value = NA_real_,
            trend_interpretation = "Model did not estimate year trend"
          )
        )
      }

      coef_year <- summary(fit)$coefficients["year_numeric", ]
      estimate <- coef_year[["Estimate"]]
      se <- coef_year[["Std. Error"]]
      p_value <- coef_year[["Pr(>|z|)"]]

      odds_ratio <- exp(estimate)
      ci_low <- exp(estimate - 1.96 * se)
      ci_high <- exp(estimate + 1.96 * se)

      tibble(
        outcome = outcome_label,
        n = nrow(dd),
        n_years = n_distinct(dd$year_numeric),
        odds_ratio_per_year = odds_ratio,
        ci_low = ci_low,
        ci_high = ci_high,
        p_value = p_value,
        trend_interpretation = case_when(
          ci_high < 1 ~ "Statistically clearer improvement",
          ci_low > 1 ~ "Statistically clearer worsening",
          odds_ratio < 1 ~ "Directionally less negative, not statistically clear",
          odds_ratio > 1 ~ "Directionally more negative, not statistically clear",
          TRUE ~ "No change"
        )
      )
    }) %>%
    ungroup()
}

sl_trend_models <- bind_rows(
  fit_municipality_trend(hhs, "food_worry", "Food worry"),
  fit_municipality_trend(hhs, "financial_strain", "Financial strain")
)

sl_trend_models %>%
  mutate(
    odds_ratio_per_year = round(odds_ratio_per_year, 3),
    ci_low = round(ci_low, 3),
    ci_high = round(ci_high, 3),
    p_value = pvalue(p_value)
  ) %>%
  arrange(outcome, odds_ratio_per_year) %>%
  kable(caption = "Exploratory municipality-specific time trends in negative Sustainable Livelihoods indicators")
Exploratory municipality-specific time trends in negative Sustainable Livelihoods indicators
g1_municipality outcome n n_years odds_ratio_per_year ci_low ci_high p_value trend_interpretation
Dondo Financial strain 213 2 0.348 0.199 0.610 <0.001 Statistically clearer improvement
Nacala Porto Financial strain 441 2 0.414 0.334 0.513 <0.001 Statistically clearer improvement
Mogincual Financial strain 614 2 0.560 0.402 0.780 <0.001 Statistically clearer improvement
Massinga Financial strain 270 2 0.840 0.658 1.072 0.161 Directionally less negative, not statistically clear
Matutuíne Financial strain 302 2 0.937 0.746 1.175 0.571 Directionally less negative, not statistically clear
Inharrime Financial strain 423 2 0.966 0.796 1.171 0.722 Directionally less negative, not statistically clear
Inhassoro Financial strain 1700 3 0.970 0.919 1.024 0.274 Directionally less negative, not statistically clear
Memba Financial strain 1302 4 1.125 1.064 1.190 <0.001 Statistically clearer worsening
Ilha de Mocambique Financial strain 1908 4 1.188 1.135 1.244 <0.001 Statistically clearer worsening
Inhassoro Food worry 1709 3 0.932 0.858 1.013 0.096 Directionally less negative, not statistically clear
Inharrime Food worry 410 2 0.983 0.777 1.243 0.884 Directionally less negative, not statistically clear
Nacala Porto Food worry 441 2 1.249 0.963 1.620 0.094 Directionally more negative, not statistically clear
Memba Food worry 1309 4 1.263 1.164 1.372 <0.001 Statistically clearer worsening
Ilha de Mocambique Food worry 1915 4 1.670 1.459 1.911 <0.001 Statistically clearer worsening
Massinga Food worry 273 2 2.434 1.656 3.578 <0.001 Statistically clearer worsening
Matutuíne Food worry 309 2 2.911 2.254 3.758 <0.001 Statistically clearer worsening
Mogincual Food worry 614 2 89.580 21.882 366.716 <0.001 Statistically clearer worsening
Dondo Food worry 226 2 NA NA NA NA Insufficient variation/data 
sl_trend_plot <- sl_trend_models %>%
  filter(!is.na(odds_ratio_per_year)) %>%
  mutate(
    municipality_plot = fct_reorder(g1_municipality, odds_ratio_per_year, .fun = stats::median, na.rm = TRUE)
  )

ggplot(
  sl_trend_plot,
  aes(
    x = odds_ratio_per_year,
    y = municipality_plot
  )
) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey40") +
  geom_segment(
    aes(x = ci_low, xend = ci_high, y = municipality_plot, yend = municipality_plot),
    alpha = 0.55
  ) +
  geom_point(aes(size = n, color = trend_interpretation), alpha = 0.85) +
  facet_wrap(~ outcome, ncol = 1, scales = "free_y") +
  scale_x_log10() +
  labs(
    title = "Municipality-specific time trends in negative livelihood conditions",
    subtitle = "Odds ratios below 1 suggest the condition is becoming less common over time; CIs crossing 1 are not statistically clear",
    x = "Odds ratio per calendar year, log scale",
    y = NULL,
    color = "Interpretation",
    size = "N"
  ) +
  theme(legend.position = "bottom")

### Interpretation note

The plot suggests that financial strain is improving in some municipalities, but food worry is not clearly improving and may be worsening in several places.

  • For financial strain, the strongest evidence of improvement is in Mogincual and Nacala Porto, where the odds ratios are below 1 and the confidence intervals appear fully below 1. This means that, in those municipalities, the odds of households reporting difficulty meeting needs have declined over time. Dondo also points toward improvement, but the confidence interval is wide, so it is not statistically clear. Ilha de Moçambique and Memba point in the opposite direction, with odds ratios above 1, suggesting financial strain may be becoming more common over time there.

  • For food worry, the pattern is less encouraging. Several municipalities have odds ratios above 1, meaning food worry appears to be becoming more common over time. This is especially clear for Ilha de Moçambique, Memba, Matutuíne, Massinga, and possibly Mogincual, although Mogincual has a very wide confidence interval, suggesting high uncertainty. Inhassoro is the main municipality that points toward improvement in food worry, but the estimate is not clearly significant.

9. Catch Data, CPUE, and artisanal-fishing income at the municipality level

This section adds a cross-data-stream exploratory analysis combining Mozambique Catch Data and HHS income-source data. The goal is to check whether municipalities and years with higher catch per unit effort (CPUE, kg/trip) also show higher HHS-derived average monthly household income from artisanal fishing.

## [1] "/Users/marianoviz/Desktop/R Projects and Stuff/ff_hhs_data_processing/data/raw/cpue_kg_trip (3).csv"
## [1] "/Users/marianoviz/Desktop/R Projects and Stuff/ff_hhs_data_processing/data/raw/join_footprint_ma.csv"
Catch Data records by country before filtering
country records
Indonesia 22886
Philippines 10849
Mozambique 7782
Honduras 2648
Guatemala 188
Managed-access mapping records by country before filtering
country_name records
Philippines 1846
Brazil 731
Indonesia 713
Honduras 92
Federated States of Micronesia 36
Mozambique 36
Guatemala 22
Palau 3 
standardize_place <- function(x) {
  x %>%
    as.character() %>%
    str_squish() %>%
    str_to_lower() %>%
    iconv(from = "UTF-8", to = "ASCII//TRANSLIT") %>%
    str_replace_all("[^a-z0-9]+", " ") %>%
    str_squish()
}

weighted_mean_ci <- function(x, w) {
  valid <- !is.na(x) & !is.na(w) & w > 0
  x <- x[valid]
  w <- w[valid]

  if (length(x) == 0) {
    return(tibble(
      mean = NA_real_,
      se = NA_real_,
      ci_low = NA_real_,
      ci_high = NA_real_,
      n_eff = NA_real_
    ))
  }

  weighted_mean <- weighted.mean(x, w, na.rm = TRUE)

  if (length(x) < 2) {
    return(tibble(
      mean = weighted_mean,
      se = NA_real_,
      ci_low = NA_real_,
      ci_high = NA_real_,
      n_eff = 1
    ))
  }

  # Approximate design-free weighted SE.
  # This treats fisher/month records as independent and uses Kish effective n.
  n_eff <- sum(w)^2 / sum(w^2)
  weighted_var <- sum(w * (x - weighted_mean)^2, na.rm = TRUE) / sum(w, na.rm = TRUE)
  se <- sqrt(weighted_var / n_eff)

  tibble(
    mean = weighted_mean,
    se = se,
    ci_low = pmax(0, weighted_mean - 1.96 * se),
    ci_high = weighted_mean + 1.96 * se,
    n_eff = n_eff
  )
}
Mozambique managed-access areas mapped to municipalities and communities
ma_key ma_name_map catch_municipality n_mapped_communities mapped_communities
gelo Gelo Angoche 1 Gelo
sangage Sangage Angoche 1 Sangage
farol Farol Dondo 1 Farol
sengo Sengo Dondo 1 Sengo
ilha insular Ilha Insular Ilha de Mocambique 1 Ilha Insular
quissanga Quissanga Ilha de Mocambique 1 Quissanga
sanculo Sanculo Ilha de Mocambique 1 Sanculo
zavora Zavora Inharrime 1 Zavora
fequete Fequete Inhassoro 1 Fequete
mucocuene Mucocuene Inhassoro 1 Mucocuene
nhagondzo Nhagondzo Inhassoro 1 Nhagondzo
petane Petane Inhassoro 1 Petane
tsondzo Tsondzo Inhassoro 1 Vuca
vuca Vuca Inhassoro 1 Tsondzo
larde sede Larde-sede Larde 1 Larde-sede
tibane Tibane Larde 1 Tibane
pomene Pomene Massinga 1 Pomene
machangulo Machangulo Matutuíne 2 Mabuluku, Santa Maria
baixo pinda Baixo Pinda Memba 1 Baixo Pinda
memba sede Memba-sede Memba 1 Memba-sede
serissa Serissa Memba 1 Serissa
simuco Simuco Memba 1 Simuco
meculuvelane Meculuvelane Mogincual 1 Meculuvelane
namalungo Namalungo Mogincual 1 Namalungo
namige sede Namige-sede Mogincual 1 Namige Sede
moma sede Moma-sede Moma 1 Moma-sede
mucoroge Mucoroge Moma 1 Mucoroge
mahelene Mahelene Nacala Porto 1 Mahelene
naherengue Naherengue Nacala Porto 1 Naherengue
quissimajulo Quissimajulo Nacala Porto 1 Quissimajulo
malaua Malaua Pebane 1 Malaua
maverane Maverane Pebane 1 Maverane
guitine Guitine Vilankulo 1 Guitine
mabandene Mabandene Vilankulo 1 Mabandene
macunhe Macunhe Vilankulo 1 Macunhe 
catch_prepared <- catch_raw %>%
  filter(str_to_lower(str_squish(country)) == "mozambique") %>%
  mutate(
    catch_year = as.integer(parse_number(as.character(year))),
    catch_month = as.integer(parse_number(as.character(month))),
    ma_name = as.character(ma_name),
    ma_key = standardize_place(ma_name),
    fisher_id = as.character(fisher_id),
    total_weight_kg = parse_number(as.character(total_weight)),
    total_price = parse_number(as.character(total_price)),
    num_trips = parse_number(as.character(num_trips)),
    cpue_raw = parse_number(as.character(cpue_kg_per_trip))
  ) %>%
  filter(
    !is.na(ma_name),
    !is.na(catch_year),
    !is.na(cpue_raw),
    cpue_raw >= 0,
    !is.na(num_trips),
    num_trips > 0
  ) %>%
  left_join(ma_crosswalk_moz, by = "ma_key")

catch_unmatched_ma <- catch_prepared %>%
  filter(is.na(catch_municipality)) %>%
  distinct(ma_name, ma_key) %>%
  arrange(ma_name)

if (nrow(catch_unmatched_ma) > 0) {
  catch_unmatched_ma %>%
    kable(caption = "Catch Data managed-access names not matched to the mapping file")
} else {
  cat("All Mozambique Catch Data managed-access names matched to the mapping file.\n")
}
## All Mozambique Catch Data managed-access names matched to the mapping file.
catch_coverage_overall <- catch_prepared %>%
  summarise(
    records = n(),
    total_trips = sum(num_trips, na.rm = TRUE),
    managed_access_areas = n_distinct(ma_name),
    mapped_municipalities = n_distinct(catch_municipality, na.rm = TRUE),
    years = paste(sort(unique(catch_year)), collapse = ", "),
    median_cpue = median(cpue_raw, na.rm = TRUE),
    mean_cpue = mean(cpue_raw, na.rm = TRUE),
    max_cpue = max(cpue_raw, na.rm = TRUE)
  )

catch_coverage_overall %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  kable(caption = "Mozambique Catch Data coverage after filtering and mapping")
Mozambique Catch Data coverage after filtering and mapping
records total_trips managed_access_areas mapped_municipalities years median_cpue mean_cpue max_cpue
7782 13714 23 10 2019, 2020, 2021, 2022, 2023, 2024, 2025, 2026 20 46.82 2000

Clean extreme CPUE values

The Catch Data file already contains cpue_kg_per_trip and num_trips. For aggregation, this analysis treats CPUE as kg/trip and weights the annual mean CPUE by num_trips. This gives more influence to records representing more trips.

Extreme CPUE values are excluded above the 99th percentile within each managed-access-year cell when the cell has at least 20 records. This cleaning is intended to reduce visual and summary distortion from extreme values while retaining small cells.

cpue_cutoffs <- catch_prepared %>%
  group_by(ma_name, catch_year) %>%
  summarise(
    n_records_cell = n(),
    total_trips_cell = sum(num_trips, na.rm = TRUE),
    cpue_p99 = if_else(
      n_records_cell >= 20,
      quantile(cpue_raw, 0.99, na.rm = TRUE),
      Inf
    ),
    .groups = "drop"
  )

catch_clean <- catch_prepared %>%
  left_join(cpue_cutoffs, by = c("ma_name", "catch_year")) %>%
  mutate(cpue_outlier_p99 = cpue_raw > cpue_p99) %>%
  filter(!cpue_outlier_p99) %>%
  mutate(cpue = cpue_raw)

catch_outlier_summary <- catch_prepared %>%
  left_join(cpue_cutoffs, by = c("ma_name", "catch_year")) %>%
  mutate(cpue_outlier_p99 = cpue_raw > cpue_p99) %>%
  group_by(catch_municipality, ma_name, catch_year) %>%
  summarise(
    records_before_cleaning = n(),
    trips_before_cleaning = sum(num_trips, na.rm = TRUE),
    records_removed_as_cpue_outliers = sum(cpue_outlier_p99, na.rm = TRUE),
    trips_removed_as_cpue_outliers = sum(num_trips[cpue_outlier_p99], na.rm = TRUE),
    .groups = "drop"
  )

catch_outlier_summary %>%
  arrange(catch_municipality, ma_name, catch_year) %>%
  kable(caption = "CPUE outlier cleaning summary by municipality, managed access area, and year")
CPUE outlier cleaning summary by municipality, managed access area, and year
catch_municipality ma_name catch_year records_before_cleaning trips_before_cleaning records_removed_as_cpue_outliers trips_removed_as_cpue_outliers
Angoche Gelo 2024 1 1 0 0
Angoche Sangage 2024 4 4 0 0
Ilha de Mocambique Ilha Insular 2020 10 24 0 0
Ilha de Mocambique Ilha Insular 2021 18 51 0 0
Ilha de Mocambique Ilha Insular 2022 402 689 4 4
Ilha de Mocambique Ilha Insular 2023 400 902 4 6
Ilha de Mocambique Ilha Insular 2024 151 332 2 3
Ilha de Mocambique Ilha Insular 2025 90 136 1 1
Ilha de Mocambique Ilha Insular 2026 18 19 0 0
Ilha de Mocambique Quissanga 2020 5 5 0 0
Ilha de Mocambique Quissanga 2021 1 1 0 0
Ilha de Mocambique Quissanga 2022 177 190 2 2
Ilha de Mocambique Quissanga 2023 82 82 1 1
Ilha de Mocambique Quissanga 2024 3 7 0 0
Ilha de Mocambique Quissanga 2025 7 10 0 0
Ilha de Mocambique Sanculo 2020 9 9 0 0
Ilha de Mocambique Sanculo 2021 6 7 0 0
Ilha de Mocambique Sanculo 2022 79 106 1 1
Ilha de Mocambique Sanculo 2023 77 95 1 1
Ilha de Mocambique Sanculo 2024 42 74 1 1
Ilha de Mocambique Sanculo 2025 32 116 1 3
Ilha de Mocambique Sanculo 2026 3 5 0 0
Inharrime Zavora 2019 59 68 1 1
Inharrime Zavora 2020 7 7 0 0
Inharrime Zavora 2022 179 482 2 7
Inharrime Zavora 2023 47 60 1 1
Inharrime Zavora 2024 2 2 0 0
Inhassoro Fequete 2019 248 383 3 7
Inhassoro Fequete 2020 193 267 2 2
Inhassoro Fequete 2021 207 287 3 6
Inhassoro Fequete 2022 310 575 4 4
Inhassoro Fequete 2023 119 203 2 5
Inhassoro Fequete 2024 100 160 1 1
Inhassoro Fequete 2025 308 438 3 3
Inhassoro Fequete 2026 361 602 4 4
Inhassoro Mucocuene 2023 4 4 0 0
Inhassoro Mucocuene 2024 11 20 0 0
Inhassoro Mucocuene 2025 38 53 1 1
Inhassoro Mucocuene 2026 53 77 1 1
Inhassoro Nhagondzo 2023 17 32 0 0
Inhassoro Nhagondzo 2024 105 127 2 2
Inhassoro Nhagondzo 2025 46 63 1 3
Inhassoro Nhagondzo 2026 51 70 1 2
Inhassoro Petane 2023 7 20 0 0
Inhassoro Petane 2024 1 1 0 0
Inhassoro Petane 2025 16 16 0 0
Inhassoro Petane 2026 11 14 0 0
Inhassoro Tsondzo 2023 52 135 1 4
Inhassoro Tsondzo 2024 65 96 1 1
Inhassoro Tsondzo 2025 27 37 1 1
Inhassoro Tsondzo 2026 11 12 0 0
Inhassoro Vuca 2023 35 54 1 2
Inhassoro Vuca 2024 5 5 0 0
Inhassoro Vuca 2025 22 24 1 1
Inhassoro Vuca 2026 11 12 0 0
Larde Larde-sede 2024 16 32 0 0
Larde Tibane 2024 2 2 0 0
Larde Tibane 2025 1 1 0 0
Massinga Pomene 2019 190 640 2 4
Massinga Pomene 2020 18 21 0 0
Massinga Pomene 2021 27 33 1 3
Massinga Pomene 2022 133 240 2 2
Massinga Pomene 2023 211 243 3 3
Massinga Pomene 2024 21 21 1 1
Matutuíne Machangulo 2019 26 28 1 1
Matutuíne Machangulo 2020 15 39 0 0
Matutuíne Machangulo 2021 1 2 0 0
Matutuíne Machangulo 2022 2 2 0 0
Matutuíne Machangulo 2023 8 11 0 0
Memba Baixo Pinda 2024 614 1617 7 22
Memba Baixo Pinda 2025 381 677 4 4
Memba Baixo Pinda 2026 18 19 0 0
Memba Memba-sede 2019 22 28 1 1
Memba Memba-sede 2020 9 9 0 0
Memba Memba-sede 2021 76 123 1 1
Memba Memba-sede 2022 134 201 2 3
Memba Memba-sede 2023 187 335 2 2
Memba Memba-sede 2024 166 257 2 2
Memba Memba-sede 2025 183 389 2 2
Memba Memba-sede 2026 33 71 1 1
Mogincual Meculuvelane 2024 54 113 1 1
Mogincual Meculuvelane 2025 33 62 1 2
Mogincual Meculuvelane 2026 8 13 0 0
Mogincual Namalungo 2024 230 342 3 3
Mogincual Namalungo 2025 71 120 1 1
Mogincual Namalungo 2026 29 39 1 2
Mogincual Namige-sede 2024 94 124 1 1
Mogincual Namige-sede 2025 78 110 1 1
Mogincual Namige-sede 2026 13 16 0 0
Nacala Porto Naherengue 2024 95 165 0 0
Nacala Porto Naherengue 2025 108 159 2 2
Nacala Porto Naherengue 2026 11 11 0 0
Nacala Porto Quissimajulo 2024 60 64 1 1
Nacala Porto Quissimajulo 2025 51 56 0 0
Nacala Porto Quissimajulo 2026 8 8 0 0

Catch-data coverage and CPUE precision by municipality-year

Precision is summarized using the standard error and approximate 95% confidence interval of the trip-weighted mean CPUE within each municipality and year. The relative margin of error is calculated as the CI half-width divided by mean CPUE. This is a descriptive precision diagnostic, not a full design-based variance estimate.

catch_precision_municipality_year <- catch_clean %>%
  filter(!is.na(catch_municipality)) %>%
  group_by(catch_municipality, catch_year) %>%
  summarise(
    n_records = n(),
    n_trips = sum(num_trips, na.rm = TRUE),
    n_unique_fishers = n_distinct(fisher_id, na.rm = TRUE),
    n_ma = n_distinct(ma_name, na.rm = TRUE),
    n_months = n_distinct(catch_month, na.rm = TRUE),
    total_weight_kg = sum(total_weight_kg, na.rm = TRUE),
    median_cpue = median(cpue, na.rm = TRUE),
    cpue_ci = list(weighted_mean_ci(cpue, num_trips)),
    .groups = "drop"
  ) %>%
  unnest(cpue_ci) %>%
  rename(
    mean_cpue = mean,
    se_cpue = se,
    ci_low_cpue = ci_low,
    ci_high_cpue = ci_high
  ) %>%
  mutate(
    relative_moe_pct = if_else(
      mean_cpue > 0 & !is.na(se_cpue),
      100 * (1.96 * se_cpue) / mean_cpue,
      NA_real_
    ),
    precision_category = case_when(
      is.na(relative_moe_pct) ~ "Not estimable",
      relative_moe_pct <= 25 ~ "Higher precision",
      relative_moe_pct <= 50 ~ "Moderate precision",
      relative_moe_pct <= 100 ~ "Low precision",
      TRUE ~ "Very low precision"
    ),
    municipality_key = standardize_place(catch_municipality)
  )

catch_precision_municipality_year %>%
  mutate(
    n_trips = round(n_trips, 0),
    total_weight_kg = round(total_weight_kg, 1),
    mean_cpue = round(mean_cpue, 2),
    median_cpue = round(median_cpue, 2),
    ci_low_cpue = round(ci_low_cpue, 2),
    ci_high_cpue = round(ci_high_cpue, 2),
    relative_moe_pct = round(relative_moe_pct, 1)
  ) %>%
  arrange(catch_municipality, catch_year) %>%
  kable(caption = "CPUE estimates and precision by municipality and year")
CPUE estimates and precision by municipality and year
catch_municipality catch_year n_records n_trips n_unique_fishers n_ma n_months total_weight_kg median_cpue mean_cpue se_cpue ci_low_cpue ci_high_cpue n_eff relative_moe_pct precision_category municipality_key
Angoche 2024 5 5 4 2 2 205.0 30.00 41.00 14.170392 13.23 68.77 5.000000 67.7 Low precision angoche
Ilha de Mocambique 2020 24 38 22 3 2 2709.5 42.00 71.30 14.452673 42.98 99.63 6.504505 39.7 Moderate precision ilha de mocambique
Ilha de Mocambique 2021 25 59 19 3 7 5268.2 20.00 89.29 23.429531 43.37 135.21 6.227191 51.4 Low precision ilha de mocambique
Ilha de Mocambique 2022 651 978 437 3 7 28790.8 20.00 29.44 1.671435 26.16 32.71 316.506949 11.1 Higher precision ilha de mocambique
Ilha de Mocambique 2023 553 1071 293 3 12 30207.8 19.50 28.21 1.116051 26.02 30.39 229.270638 7.8 Higher precision ilha de mocambique
Ilha de Mocambique 2024 193 409 113 3 11 11548.3 17.00 28.24 2.599331 23.14 33.33 83.515227 18.0 Higher precision ilha de mocambique
Ilha de Mocambique 2025 127 258 56 3 12 4989.3 13.83 19.34 1.846627 15.72 22.96 61.405904 18.7 Higher precision ilha de mocambique
Ilha de Mocambique 2026 21 24 17 2 5 225.0 3.00 9.38 3.905409 1.72 17.03 18.000000 81.6 Low precision ilha de mocambique
Inharrime 2019 58 67 34 1 6 2101.5 19.00 31.37 5.036252 21.49 41.24 45.343434 31.5 Moderate precision inharrime
Inharrime 2020 7 7 7 1 1 629.0 75.00 89.86 17.426731 55.70 124.01 7.000000 38.0 Moderate precision inharrime
Inharrime 2022 177 475 68 1 4 11166.0 18.50 23.51 1.162315 21.23 25.79 112.868934 9.7 Higher precision inharrime
Inharrime 2023 46 59 27 1 8 1823.5 21.25 30.91 3.482206 24.08 37.73 37.430107 22.1 Higher precision inharrime
Inharrime 2024 2 2 2 1 2 16.0 8.00 8.00 1.414214 5.23 10.77 2.000000 34.6 Moderate precision inharrime
Inhassoro 2019 245 376 139 1 10 17113.0 20.00 45.51 4.166682 37.35 53.68 170.743961 17.9 Higher precision inhassoro
Inhassoro 2020 191 265 131 1 12 13005.7 25.00 49.08 5.723320 37.86 60.30 120.042735 22.9 Higher precision inhassoro
Inhassoro 2021 204 281 130 1 12 24114.3 35.62 85.82 8.537971 69.08 102.55 127.975689 19.5 Higher precision inhassoro
Inhassoro 2022 306 571 150 1 12 43221.8 35.00 75.69 6.639419 62.68 88.71 157.583857 17.2 Higher precision inhassoro
Inhassoro 2023 230 437 163 6 7 13956.2 10.70 31.94 5.059561 22.02 41.85 138.083153 31.1 Moderate precision inhassoro
Inhassoro 2024 283 405 147 6 12 10022.0 12.00 24.75 2.674723 19.50 29.99 200.764994 21.2 Higher precision inhassoro
Inhassoro 2025 450 622 256 6 12 65885.0 50.00 105.92 8.137665 89.97 121.87 298.982999 15.1 Higher precision inhassoro
Inhassoro 2026 492 780 286 6 7 61954.2 30.00 79.43 5.882014 67.90 90.96 276.797088 14.5 Higher precision inhassoro
Larde 2024 18 34 17 2 2 965.0 20.00 28.38 5.080370 18.42 38.34 13.761905 35.1 Moderate precision larde
Larde 2025 1 1 1 1 1 25.0 25.00 25.00 NA NA NA 1.000000 NA Not estimable larde
Massinga 2019 188 636 56 1 9 9064.0 12.00 14.25 0.949332 12.39 16.11 115.438356 13.1 Higher precision massinga
Massinga 2020 18 21 15 1 3 308.5 13.00 14.69 1.762698 11.24 18.15 15.206897 23.5 Higher precision massinga
Massinga 2021 26 30 21 1 7 515.0 15.75 17.17 1.850244 13.54 20.79 23.684210 21.1 Higher precision massinga
Massinga 2022 131 238 76 1 5 9165.5 30.00 38.51 2.509992 33.59 43.43 75.324468 12.8 Higher precision massinga
Massinga 2023 208 240 85 1 12 13368.0 39.00 55.70 3.454901 48.93 62.47 180.000000 12.2 Higher precision massinga
Massinga 2024 20 20 18 1 3 1468.0 65.50 73.40 8.621891 56.50 90.30 20.000000 23.0 Higher precision massinga
Matutuíne 2019 25 27 15 1 9 1155.0 38.00 42.78 4.587045 33.79 51.77 23.516129 21.0 Higher precision matutu ine
Matutuíne 2020 15 39 8 1 9 1796.5 30.00 46.06 11.903482 22.73 69.39 6.059761 50.6 Low precision matutu ine
Matutuíne 2021 1 2 1 1 1 277.0 138.50 138.50 NA NA NA 1.000000 NA Not estimable matutu ine
Matutuíne 2022 2 2 2 1 1 26.0 13.00 13.00 4.949747 3.30 22.70 2.000000 74.6 Low precision matutu ine
Matutuíne 2023 8 11 6 1 2 331.5 20.12 30.14 12.239783 6.15 54.13 7.117647 79.6 Low precision matutu ine
Memba 2019 21 27 18 1 5 221.0 5.00 8.19 1.352364 5.53 10.84 13.754717 32.4 Moderate precision memba
Memba 2020 9 9 9 1 4 750.0 62.00 83.33 25.628880 33.10 133.57 9.000000 60.3 Low precision memba
Memba 2021 75 122 42 1 9 1890.0 10.00 15.49 3.378284 8.87 22.11 45.938272 42.7 Moderate precision memba
Memba 2022 132 198 77 1 9 6283.3 19.00 31.73 3.688452 24.50 38.96 83.059322 22.8 Higher precision memba
Memba 2023 185 333 56 1 12 11409.0 20.00 34.26 4.007595 26.41 42.12 101.826446 22.9 Higher precision memba
Memba 2024 771 1850 263 2 12 39547.3 13.00 21.38 1.029880 19.36 23.40 409.389952 9.4 Higher precision memba
Memba 2025 558 1060 218 2 12 20025.4 12.00 18.89 1.016809 16.90 20.88 315.263749 10.5 Higher precision memba
Memba 2026 50 89 37 2 5 1361.0 10.00 15.29 2.220345 10.94 19.64 31.811245 28.5 Moderate precision memba
Mogincual 2024 373 574 218 3 7 58306.0 60.00 101.58 9.294226 83.36 119.80 225.359781 17.9 Higher precision mogincual
Mogincual 2025 179 288 68 3 12 16047.0 43.00 55.72 4.117312 47.65 63.79 120.208696 14.5 Higher precision mogincual
Mogincual 2026 49 66 30 3 6 3567.0 50.00 54.05 6.413076 41.48 66.62 41.094340 23.3 Higher precision mogincual
Nacala Porto 2024 154 228 89 2 12 7493.5 18.79 32.87 4.481942 24.08 41.65 99.206107 26.7 Moderate precision nacala porto
Nacala Porto 2025 157 213 96 2 12 7685.5 25.00 36.08 5.728358 24.85 47.31 110.926650 31.1 Moderate precision nacala porto
Nacala Porto 2026 19 19 16 2 7 1154.0 40.00 60.74 16.481811 28.43 93.04 19.000000 53.2 Low precision nacala porto 
catch_trip_heatmap <- catch_precision_municipality_year %>%
  mutate(
    municipality_plot = fct_reorder(catch_municipality, n_trips, .fun = sum),
    label = paste0("trips=", comma(round(n_trips, 0)), "\nrecords=", n_records),
    value_scaled = rescale(n_trips, to = c(0, 1), from = range(n_trips, na.rm = TRUE)),
    label_color = if_else(value_scaled < 0.45, "white", "black")
  )

catch_trip_heatmap %>%
  ggplot(
    aes(
      x = factor(catch_year),
      y = municipality_plot,
      fill = n_trips
    )
  ) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 2.8, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(labels = comma, option = "C") +
  labs(
    title = "Mozambique Catch Data coverage by municipality and year",
    subtitle = "Cell labels show valid trips and records after CPUE outlier cleaning",
    x = "Year",
    y = NULL,
    fill = "Trips"
  ) +
  theme(panel.grid = element_blank())

precision_heatmap <- catch_precision_municipality_year %>%
  mutate(
    municipality_plot = fct_reorder(catch_municipality, relative_moe_pct, .fun = stats::median, na.rm = TRUE),
    label = case_when(
      is.na(relative_moe_pct) ~ paste0("n=", n_records, "\nNA"),
      TRUE ~ paste0("RMoE=", round(relative_moe_pct, 0), "%", "\ntrips=", comma(round(n_trips, 0)))
    ),
    value_scaled = rescale(relative_moe_pct, to = c(0, 1), from = range(relative_moe_pct, na.rm = TRUE)),
    label_color = if_else(value_scaled > 0.60, "white", "black")
  )

precision_heatmap %>%
  ggplot(
    aes(
      x = factor(catch_year),
      y = municipality_plot,
      fill = relative_moe_pct
    )
  ) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 2.8, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(
    labels = label_percent(scale = 1),
    option = "C",
    direction = -1,
    na.value = "grey90"
  ) +
  labs(
    title = "Statistical precision of CPUE estimates by municipality and year",
    subtitle = "Relative margin of error = 95% CI half-width divided by trip-weighted mean CPUE; lower is better",
    x = "Year",
    y = NULL,
    fill = "Relative MOE"
  ) +
  theme(panel.grid = element_blank())

ggplot(
  catch_precision_municipality_year,
  aes(
    x = catch_year,
    y = mean_cpue,
    group = catch_municipality,
    color = catch_municipality
  )
) +
  geom_errorbar(aes(ymin = ci_low_cpue, ymax = ci_high_cpue), width = 0.12, alpha = 0.45) +
  geom_line(linewidth = 0.8) +
  geom_point(aes(size = n_trips), alpha = 0.85) +
  scale_x_continuous(breaks = sort(unique(catch_precision_municipality_year$catch_year))) +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Trip-weighted mean CPUE by municipality and year",
    subtitle = "Error bars show approximate 95% CIs; point size reflects number of trips",
    x = "Year",
    y = "Mean CPUE, kg/trip",
    color = "Municipality",
    size = "Trips"
  ) +
  theme(legend.position = "bottom")

HHS-derived average household income from artisanal fishing by municipality-year

For the income side, the analysis calculates average monthly household income from artisanal fishing as:

\[ \text{artisanal fishing income}_{2021\ MZN} = \text{real monthly household income}_{2021\ MZN} \times \frac{\text{% income from artisanal fishing}}{100} \]

The income variable is the already cleaned real household income variable, where values above the year-specific p99 are excluded.

hhs_fishing_income_municipality_year <- hhs %>%
  mutate(
    municipality_key = standardize_place(g1_municipality),
    hh_artisanal_fishing_income_2021_mzn =
      real_monthly_hh_income_2021_mzn_clean * (fishing_artisanal_income_pct / 100)
  ) %>%
  filter(
    !is.na(year),
    !is.na(g1_municipality),
    !is.na(hh_artisanal_fishing_income_2021_mzn),
    hh_artisanal_fishing_income_2021_mzn >= 0
  ) %>%
  group_by(g1_municipality, municipality_key, year) %>%
  summarise(
    n_hhs = n(),
    mean_hh_artisanal_fishing_income = mean(hh_artisanal_fishing_income_2021_mzn, na.rm = TRUE),
    sd_hh_artisanal_fishing_income = sd(hh_artisanal_fishing_income_2021_mzn, na.rm = TRUE),
    se_hh_artisanal_fishing_income = sd_hh_artisanal_fishing_income / sqrt(n_hhs),
    ci_low_hh_artisanal_fishing_income = pmax(0, mean_hh_artisanal_fishing_income - 1.96 * se_hh_artisanal_fishing_income),
    ci_high_hh_artisanal_fishing_income = mean_hh_artisanal_fishing_income + 1.96 * se_hh_artisanal_fishing_income,
    median_hh_artisanal_fishing_income = median(hh_artisanal_fishing_income_2021_mzn, na.rm = TRUE),
    q25_hh_artisanal_fishing_income = quantile(hh_artisanal_fishing_income_2021_mzn, 0.25, na.rm = TRUE),
    q75_hh_artisanal_fishing_income = quantile(hh_artisanal_fishing_income_2021_mzn, 0.75, na.rm = TRUE),
    mean_artisanal_fishing_dependency_pct = mean(fishing_artisanal_income_pct, na.rm = TRUE),
    .groups = "drop"
  )

hhs_fishing_income_municipality_year %>%
  mutate(
    mean_hh_artisanal_fishing_income = round(mean_hh_artisanal_fishing_income, 0),
    median_hh_artisanal_fishing_income = round(median_hh_artisanal_fishing_income, 0),
    mean_artisanal_fishing_dependency_pct = round(mean_artisanal_fishing_dependency_pct, 1)
  ) %>%
  arrange(g1_municipality, year) %>%
  kable(caption = "HHS-derived average monthly household income from artisanal fishing by municipality and year")
HHS-derived average monthly household income from artisanal fishing by municipality and year
g1_municipality municipality_key year n_hhs mean_hh_artisanal_fishing_income sd_hh_artisanal_fishing_income se_hh_artisanal_fishing_income ci_low_hh_artisanal_fishing_income ci_high_hh_artisanal_fishing_income median_hh_artisanal_fishing_income q25_hh_artisanal_fishing_income q75_hh_artisanal_fishing_income mean_artisanal_fishing_dependency_pct
Dondo dondo 2019 189 0 0.0000 0.00000 0.0000 0.0000 0 0.000 0.000 0.0
Dondo dondo 2021 9 0 0.0000 0.00000 0.0000 0.0000 0 0.000 0.000 0.0
Ilha de Mocambique ilha de mocambique 2019 325 2220 4055.2065 224.94239 1779.4235 2661.1977 0 0.000 3523.840 31.0
Ilha de Mocambique ilha de mocambique 2021 298 2322 2932.9601 169.90182 1989.0847 2655.0998 1000 0.000 4000.000 61.6
Ilha de Mocambique ilha de mocambique 2023 309 5177 4990.1927 283.88214 4620.4083 5733.2263 4063 1523.700 7237.575 33.5
Ilha de Mocambique ilha de mocambique 2025 941 2950 5254.4381 171.28979 2614.4639 3285.9199 1091 0.000 3584.320 43.2
Inharrime inharrime 2019 186 2615 4495.5789 329.63184 1968.5315 3260.6884 0 0.000 4404.800 27.9
Inharrime inharrime 2021 206 4380 5676.4660 395.49813 3604.7023 5155.0550 0 0.000 8537.500 37.7
Inhassoro inhassoro 2019 150 7065 8636.6088 705.17616 5682.8236 8447.1141 2588 0.000 13214.400 50.9
Inhassoro inhassoro 2021 843 1643 3193.6189 109.99408 1427.3950 1858.5718 0 0.000 2450.000 30.1
Inhassoro inhassoro 2025 598 2202 3378.9234 138.17445 1931.1269 2472.7708 573 0.000 3740.160 31.6
Massinga massinga 2019 81 8659 8790.1806 976.68674 6744.4160 10573.0280 4955 110.120 16518.000 70.1
Massinga massinga 2021 131 4216 6346.6710 554.51122 3129.4176 5303.1015 0 0.000 6000.000 39.5
Matutuíne matutu ine 2019 144 5370 4757.3807 396.44839 4592.9935 6147.0712 4405 1527.915 8011.230 58.3
Matutuíne matutu ine 2021 135 1813 1430.1297 123.08596 1572.0348 2054.5318 2125 475.000 2500.000 24.5
Memba memba 2019 154 771 1386.1302 111.69756 552.1530 990.0075 0 0.000 991.080 25.5
Memba memba 2021 652 225 505.1768 19.78425 186.1919 263.7462 0 0.000 360.000 27.5
Memba memba 2024 145 1822 919.2382 76.33858 1672.8134 1972.0606 2033 1219.950 2439.900 54.1
Memba memba 2026 274 616 1084.6783 65.52781 487.1513 744.0203 358 104.496 783.720 38.4
Mogincual mogincual 2024 298 10836 27055.9971 1567.31183 7763.6715 13907.5339 2440 0.000 7319.700 39.5
Mogincual mogincual 2025 308 1738 2873.3721 163.72557 1416.7348 2058.5390 1122 496.740 2175.916 41.5
Nacala Porto nacala porto 2024 260 2611 4779.8207 296.43190 2030.2489 3192.2619 651 0.000 3304.031 45.8
Nacala Porto nacala porto 2026 177 3862 5556.5736 417.65764 3043.5020 4680.7200 1493 0.000 5821.920 40.7 
hhs_fishing_income_heatmap <- hhs_fishing_income_municipality_year %>%
  mutate(
    municipality_plot = fct_reorder(g1_municipality, mean_hh_artisanal_fishing_income, .fun = stats::median, na.rm = TRUE),
    label = paste0(comma(round(mean_hh_artisanal_fishing_income, 0)), "\n", "n=", n_hhs),
    value_scaled = rescale(
      mean_hh_artisanal_fishing_income,
      to = c(0, 1),
      from = range(mean_hh_artisanal_fishing_income, na.rm = TRUE)
    ),
    label_color = if_else(value_scaled < 0.45, "white", "black")
  )

hhs_fishing_income_heatmap %>%
  ggplot(
    aes(
      x = factor(year),
      y = municipality_plot,
      fill = mean_hh_artisanal_fishing_income
    )
  ) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 3, fontface = "bold") +
  scale_color_identity() +
  scale_fill_viridis_c(labels = comma, option = "C") +
  labs(
    title = "Average household income from artisanal fishing by municipality and year",
    subtitle = "Income shown in constant 2021 MZN; cell labels show mean income and HHS sample size",
    x = "Survey year",
    y = NULL,
    fill = "Mean income"
  ) +
  theme(panel.grid = element_blank())

Match municipality-year CPUE to HHS-derived artisanal-fishing income

The join uses standardized municipality names from Catch Data mapping (lgu_name) and HHS (g1_municipality), plus year. Only municipality-years that have both Catch Data and HHS-derived income data can be used in the CPUE-income relationship analysis.

cpue_income_municipality_year <- catch_precision_municipality_year %>%
  left_join(
    hhs_fishing_income_municipality_year,
    by = c("municipality_key" = "municipality_key", "catch_year" = "year")
  ) %>%
  mutate(
    matched_hhs = !is.na(mean_hh_artisanal_fishing_income)
  )

matched_cpue_income <- cpue_income_municipality_year %>%
  filter(matched_hhs)

cpue_income_municipality_year %>%
  mutate(
    mean_cpue = round(mean_cpue, 2),
    ci_low_cpue = round(ci_low_cpue, 2),
    ci_high_cpue = round(ci_high_cpue, 2),
    relative_moe_pct = round(relative_moe_pct, 1),
    mean_hh_artisanal_fishing_income = round(mean_hh_artisanal_fishing_income, 0),
    median_hh_artisanal_fishing_income = round(median_hh_artisanal_fishing_income, 0)
  ) %>%
  arrange(catch_municipality, catch_year) %>%
  select(
    catch_municipality,
    catch_year,
    n_ma,
    n_trips,
    n_records,
    mean_cpue,
    ci_low_cpue,
    ci_high_cpue,
    relative_moe_pct,
    g1_municipality,
    n_hhs,
    mean_hh_artisanal_fishing_income,
    median_hh_artisanal_fishing_income,
    matched_hhs
  ) %>%
  kable(caption = "Matched CPUE and HHS-derived artisanal-fishing income by municipality and year")
Matched CPUE and HHS-derived artisanal-fishing income by municipality and year
catch_municipality catch_year n_ma n_trips n_records mean_cpue ci_low_cpue ci_high_cpue relative_moe_pct g1_municipality n_hhs mean_hh_artisanal_fishing_income median_hh_artisanal_fishing_income matched_hhs
Angoche 2024 2 5 5 41.00 13.23 68.77 67.7 NA NA NA NA FALSE
Ilha de Mocambique 2020 3 38 24 71.30 42.98 99.63 39.7 NA NA NA NA FALSE
Ilha de Mocambique 2021 3 59 25 89.29 43.37 135.21 51.4 Ilha de Mocambique 298 2322 1000 TRUE
Ilha de Mocambique 2022 3 978 651 29.44 26.16 32.71 11.1 NA NA NA NA FALSE
Ilha de Mocambique 2023 3 1071 553 28.21 26.02 30.39 7.8 Ilha de Mocambique 309 5177 4063 TRUE
Ilha de Mocambique 2024 3 409 193 28.24 23.14 33.33 18.0 NA NA NA NA FALSE
Ilha de Mocambique 2025 3 258 127 19.34 15.72 22.96 18.7 Ilha de Mocambique 941 2950 1091 TRUE
Ilha de Mocambique 2026 2 24 21 9.38 1.72 17.03 81.6 NA NA NA NA FALSE
Inharrime 2019 1 67 58 31.37 21.49 41.24 31.5 Inharrime 186 2615 0 TRUE
Inharrime 2020 1 7 7 89.86 55.70 124.01 38.0 NA NA NA NA FALSE
Inharrime 2022 1 475 177 23.51 21.23 25.79 9.7 NA NA NA NA FALSE
Inharrime 2023 1 59 46 30.91 24.08 37.73 22.1 NA NA NA NA FALSE
Inharrime 2024 1 2 2 8.00 5.23 10.77 34.6 NA NA NA NA FALSE
Inhassoro 2019 1 376 245 45.51 37.35 53.68 17.9 Inhassoro 150 7065 2588 TRUE
Inhassoro 2020 1 265 191 49.08 37.86 60.30 22.9 NA NA NA NA FALSE
Inhassoro 2021 1 281 204 85.82 69.08 102.55 19.5 Inhassoro 843 1643 0 TRUE
Inhassoro 2022 1 571 306 75.69 62.68 88.71 17.2 NA NA NA NA FALSE
Inhassoro 2023 6 437 230 31.94 22.02 41.85 31.1 NA NA NA NA FALSE
Inhassoro 2024 6 405 283 24.75 19.50 29.99 21.2 NA NA NA NA FALSE
Inhassoro 2025 6 622 450 105.92 89.97 121.87 15.1 Inhassoro 598 2202 573 TRUE
Inhassoro 2026 6 780 492 79.43 67.90 90.96 14.5 NA NA NA NA FALSE
Larde 2024 2 34 18 28.38 18.42 38.34 35.1 NA NA NA NA FALSE
Larde 2025 1 1 1 25.00 NA NA NA NA NA NA NA FALSE
Massinga 2019 1 636 188 14.25 12.39 16.11 13.1 Massinga 81 8659 4955 TRUE
Massinga 2020 1 21 18 14.69 11.24 18.15 23.5 NA NA NA NA FALSE
Massinga 2021 1 30 26 17.17 13.54 20.79 21.1 Massinga 131 4216 0 TRUE
Massinga 2022 1 238 131 38.51 33.59 43.43 12.8 NA NA NA NA FALSE
Massinga 2023 1 240 208 55.70 48.93 62.47 12.2 NA NA NA NA FALSE
Massinga 2024 1 20 20 73.40 56.50 90.30 23.0 NA NA NA NA FALSE
Matutuíne 2019 1 27 25 42.78 33.79 51.77 21.0 Matutuíne 144 5370 4405 TRUE
Matutuíne 2020 1 39 15 46.06 22.73 69.39 50.6 NA NA NA NA FALSE
Matutuíne 2021 1 2 1 138.50 NA NA NA Matutuíne 135 1813 2125 TRUE
Matutuíne 2022 1 2 2 13.00 3.30 22.70 74.6 NA NA NA NA FALSE
Matutuíne 2023 1 11 8 30.14 6.15 54.13 79.6 NA NA NA NA FALSE
Memba 2019 1 27 21 8.19 5.53 10.84 32.4 Memba 154 771 0 TRUE
Memba 2020 1 9 9 83.33 33.10 133.57 60.3 NA NA NA NA FALSE
Memba 2021 1 122 75 15.49 8.87 22.11 42.7 Memba 652 225 0 TRUE
Memba 2022 1 198 132 31.73 24.50 38.96 22.8 NA NA NA NA FALSE
Memba 2023 1 333 185 34.26 26.41 42.12 22.9 NA NA NA NA FALSE
Memba 2024 2 1850 771 21.38 19.36 23.40 9.4 Memba 145 1822 2033 TRUE
Memba 2025 2 1060 558 18.89 16.90 20.88 10.5 NA NA NA NA FALSE
Memba 2026 2 89 50 15.29 10.94 19.64 28.5 Memba 274 616 358 TRUE
Mogincual 2024 3 574 373 101.58 83.36 119.80 17.9 Mogincual 298 10836 2440 TRUE
Mogincual 2025 3 288 179 55.72 47.65 63.79 14.5 Mogincual 308 1738 1122 TRUE
Mogincual 2026 3 66 49 54.05 41.48 66.62 23.3 NA NA NA NA FALSE
Nacala Porto 2024 2 228 154 32.87 24.08 41.65 26.7 Nacala Porto 260 2611 651 TRUE
Nacala Porto 2025 2 213 157 36.08 24.85 47.31 31.1 NA NA NA NA FALSE
Nacala Porto 2026 2 19 19 60.74 28.43 93.04 53.2 Nacala Porto 177 3862 1493 TRUE 
unmatched_catch_municipality_year <- cpue_income_municipality_year %>%
  filter(!matched_hhs) %>%
  select(catch_municipality, catch_year, n_trips, n_records, mean_cpue)

if (nrow(unmatched_catch_municipality_year) > 0) {
  unmatched_catch_municipality_year %>%
    arrange(catch_municipality, catch_year) %>%
    kable(caption = "Catch municipality-years not matched to HHS income data")
}
Catch municipality-years not matched to HHS income data
catch_municipality catch_year n_trips n_records mean_cpue
Angoche 2024 5 5 41.00000
Ilha de Mocambique 2020 38 24 71.30263
Ilha de Mocambique 2022 978 651 29.43845
Ilha de Mocambique 2024 409 193 28.23545
Ilha de Mocambique 2026 24 21 9.37500
Inharrime 2020 7 7 89.85714
Inharrime 2022 475 177 23.50737
Inharrime 2023 59 46 30.90678
Inharrime 2024 2 2 8.00000
Inhassoro 2020 265 191 49.07811
Inhassoro 2022 571 306 75.69492
Inhassoro 2023 437 230 31.93638
Inhassoro 2024 405 283 24.74568
Inhassoro 2026 780 492 79.42846
Larde 2024 34 18 28.38235
Larde 2025 1 1 25.00000
Massinga 2020 21 18 14.69048
Massinga 2022 238 131 38.51050
Massinga 2023 240 208 55.70000
Massinga 2024 20 20 73.40000
Matutuíne 2020 39 15 46.06410
Matutuíne 2022 2 2 13.00000
Matutuíne 2023 11 8 30.13636
Memba 2020 9 9 83.33333
Memba 2022 198 132 31.73384
Memba 2023 333 185 34.26126
Memba 2025 1060 558 18.89189
Mogincual 2026 66 49 54.04545
Nacala Porto 2025 213 157 36.08216 
cat("Matched municipality-years available for CPUE-income analysis:", nrow(matched_cpue_income), "\n")
## Matched municipality-years available for CPUE-income analysis: 19

CPUE versus average household income from artisanal fishing

The primary relationship plot compares trip-weighted mean CPUE with mean HHS-derived monthly household income from artisanal fishing. Mean income is used here because the research question is about average household monthly fishing income, but the median is also shown as a robustness check because income is skewed.

if (nrow(matched_cpue_income) == 0) {
  message("No matched CPUE-HHS municipality-years available. Check municipality names and years.")
} else {
  ggplot(
    matched_cpue_income,
    aes(
      x = mean_cpue,
      y = mean_hh_artisanal_fishing_income
    )
  ) +
    geom_segment(
      aes(
        x = ci_low_cpue,
        xend = ci_high_cpue,
        y = mean_hh_artisanal_fishing_income,
        yend = mean_hh_artisanal_fishing_income
      ),
      alpha = 0.35
    ) +
    geom_segment(
      aes(
        x = mean_cpue,
        xend = mean_cpue,
        y = ci_low_hh_artisanal_fishing_income,
        yend = ci_high_hh_artisanal_fishing_income
      ),
      alpha = 0.35
    ) +
    geom_point(
      aes(size = pmin(n_trips, n_hhs), color = factor(catch_year)),
      alpha = 0.85
    ) +
    geom_smooth(method = "lm", se = TRUE, color = "grey30", linewidth = 0.8) +
    geom_text(
      aes(label = paste0(catch_municipality, "\n", catch_year)),
      size = 3,
      check_overlap = TRUE,
      show.legend = FALSE,
      vjust = -0.7
    ) +
    scale_x_continuous(labels = comma) +
    scale_y_continuous(labels = comma) +
    labs(
      title = "CPUE versus average household income from artisanal fishing",
      subtitle = "Each point is a matched municipality-year; horizontal bars show CPUE CI; vertical bars show HHS mean-income CI",
      x = "Trip-weighted mean CPUE, kg/trip",
      y = "Mean monthly household income from artisanal fishing, 2021 MZN",
      color = "Year",
      size = "Min(n trips, n HHS)"
    ) +
    theme(legend.position = "bottom")
}

The figure does not show a clear relationship between CPUE and average household income from artisanal fishing at the municipality-year level. The trend is slightly positive, but the points are widely scattered and uncertainty is high, suggesting that CPUE alone does not explain variation in fishing income. This may reflect limited matched data, measurement uncertainty, and the role of other factors such as prices, effort, market access, and household livelihood strategies.

if (nrow(matched_cpue_income) == 0) {
  message("No matched CPUE-HHS municipality-years available.")
} else {
  ggplot(
    matched_cpue_income,
    aes(
      x = mean_cpue,
      y = median_hh_artisanal_fishing_income
    )
  ) +
    geom_segment(
      aes(
        x = ci_low_cpue,
        xend = ci_high_cpue,
        y = median_hh_artisanal_fishing_income,
        yend = median_hh_artisanal_fishing_income
      ),
      alpha = 0.35
    ) +
    geom_segment(
      aes(
        x = mean_cpue,
        xend = mean_cpue,
        y = q25_hh_artisanal_fishing_income,
        yend = q75_hh_artisanal_fishing_income
      ),
      alpha = 0.35
    ) +
    geom_point(
      aes(size = pmin(n_trips, n_hhs), color = factor(catch_year)),
      alpha = 0.85
    ) +
    geom_smooth(method = "lm", se = TRUE, color = "grey30", linewidth = 0.8) +
    geom_text(
      aes(label = paste0(catch_municipality, "\n", catch_year)),
      size = 3,
      check_overlap = TRUE,
      show.legend = FALSE,
      vjust = -0.7
    ) +
    scale_x_continuous(labels = comma) +
    scale_y_continuous(labels = comma) +
    labs(
      title = "CPUE versus median household income from artisanal fishing",
      subtitle = "Median shown as a robustness check because income is skewed; vertical bars show IQR",
      x = "Trip-weighted mean CPUE, kg/trip",
      y = "Median monthly household income from artisanal fishing, 2021 MZN",
      color = "Year",
      size = "Min(n trips, n HHS)"
    ) +
    theme(legend.position = "bottom")
}

The median-income version reinforces the same conclusion: there is no clear relationship between CPUE and typical household income from artisanal fishing at the municipality-year level.

if (nrow(matched_cpue_income) == 0) {
  message("No matched CPUE-HHS municipality-years available.")
} else {
  ggplot(
    matched_cpue_income,
    aes(
      x = log1p(mean_cpue),
      y = log1p(mean_hh_artisanal_fishing_income)
    )
  ) +
    geom_point(aes(size = pmin(n_trips, n_hhs), color = factor(catch_year)), alpha = 0.85) +
    geom_smooth(method = "lm", se = TRUE, color = "grey30", linewidth = 0.8) +
    scale_x_continuous(labels = label_number()) +
    scale_y_continuous(labels = label_number()) +
    labs(
      title = "Log-scale relationship between CPUE and average artisanal-fishing income",
      subtitle = "The log scale reduces the influence of high-income or high-CPUE municipality-years",
      x = "log(1 + trip-weighted mean CPUE)",
      y = "log(1 + average monthly household income from artisanal fishing)",
      color = "Year",
      size = "Min(n trips, n HHS)"
    ) +
    theme(legend.position = "bottom")
}

The log-scale plot also suggests only a weak and uncertain positive relationship between CPUE and average household income from artisanal fishing. The fitted line slopes upward, meaning higher CPUE may be associated with higher fishing income, but the confidence band is wide and the municipality-year points remain quite dispersed.

Statistical tests of the CPUE-income relationship

The main test uses matched municipality-year observations. The primary outcome is log(1 + average monthly household income from artisanal fishing) and the main predictor is log(1 + trip-weighted mean CPUE). Models are fitted only if there are enough matched municipality-years with variation in CPUE.

The model set is intentionally simple because the sample size at the municipality-year level may be limited. The year-adjusted and municipality-adjusted specifications are included only when there are enough observations to support them.

fit_cpue_income_models <- function(data) {
  data <- data %>%
    filter(
      !is.na(mean_cpue),
      !is.na(mean_hh_artisanal_fishing_income),
      mean_cpue >= 0,
      mean_hh_artisanal_fishing_income >= 0
    ) %>%
    mutate(
      log_mean_cpue = log1p(mean_cpue),
      log_mean_fishing_income = log1p(mean_hh_artisanal_fishing_income),
      log_median_fishing_income = log1p(median_hh_artisanal_fishing_income),
      model_weight = sqrt(pmax(n_trips, 1) * pmax(n_hhs, 1))
    )

  if (nrow(data) < 5 || n_distinct(data$mean_cpue) < 2) {
    return(tibble(
      outcome = character(),
      model = character(),
      n_municipality_years = integer(),
      estimate_log_cpue = numeric(),
      conf_low = numeric(),
      conf_high = numeric(),
      p_value = numeric(),
      interpretation = character()
    ))
  }

  model_specs <- list(
    "Mean income | Unweighted" = lm(log_mean_fishing_income ~ log_mean_cpue, data = data),
    "Mean income | Weighted by HHS and trip coverage" = lm(log_mean_fishing_income ~ log_mean_cpue, data = data, weights = model_weight),
    "Median income | Unweighted" = lm(log_median_fishing_income ~ log_mean_cpue, data = data)
  )

  if (nrow(data) >= 8 && n_distinct(data$catch_year) > 1) {
    model_specs[["Mean income | Year-adjusted, unweighted"]] <- lm(
      log_mean_fishing_income ~ log_mean_cpue + factor(catch_year),
      data = data
    )
  }

  if (nrow(data) >= 12 && n_distinct(data$catch_municipality) > 1) {
    model_specs[["Mean income | Municipality-adjusted, unweighted"]] <- lm(
      log_mean_fishing_income ~ log_mean_cpue + factor(catch_municipality),
      data = data
    )
  }

  imap_dfr(
    model_specs,
    ~ tidy(.x, conf.int = TRUE) %>%
      filter(term == "log_mean_cpue") %>%
      transmute(
        outcome = str_extract(.y, "^[^|]+") %>% str_squish(),
        model = str_extract(.y, "(?<=\\| ).+$") %>% str_squish(),
        n_municipality_years = nobs(.x),
        estimate_log_cpue = estimate,
        conf_low = conf.low,
        conf_high = conf.high,
        p_value = p.value,
        interpretation = case_when(
          p_value < 0.05 & estimate > 0 ~ "Positive association, statistically significant at p<0.05",
          p_value < 0.05 & estimate < 0 ~ "Negative association, statistically significant at p<0.05",
          TRUE ~ "No statistically significant association at p<0.05"
        )
      )
  )
}

cpue_income_model_results <- fit_cpue_income_models(matched_cpue_income)

if (nrow(cpue_income_model_results) == 0) {
  message("Not enough matched municipality-years to fit a CPUE-income model. This is itself an important data-coverage finding.")
} else {
  cpue_income_model_results %>%
    mutate(
      estimate_log_cpue = round(estimate_log_cpue, 3),
      conf_low = round(conf_low, 3),
      conf_high = round(conf_high, 3),
      p_value = pvalue(p_value)
    ) %>%
    kable(caption = "Exploratory models: relationship between CPUE and HHS-derived artisanal-fishing income at municipality-year level")
}
Exploratory models: relationship between CPUE and HHS-derived artisanal-fishing income at municipality-year level
outcome model n_municipality_years estimate_log_cpue conf_low conf_high p_value interpretation
Mean income Unweighted 19 0.364 -0.219 0.946 0.205 No statistically significant association at p<0.05
Mean income Weighted by HHS and trip coverage 19 0.325 -0.292 0.942 0.282 No statistically significant association at p<0.05
Median income Unweighted 19 1.459 -0.614 3.532 0.156 No statistically significant association at p<0.05
Mean income Year-adjusted, unweighted 19 0.593 -0.019 1.205 0.056 No statistically significant association at p<0.05
Mean income Municipality-adjusted, unweighted 19 -0.253 -1.254 0.748 0.586 No statistically significant association at p<0.05 
if (nrow(matched_cpue_income) >= 5) {
  cpue_income_correlations <- matched_cpue_income %>%
    summarise(
      n_municipality_years = n(),
      pearson_log_cpue_mean_income = cor(
        log1p(mean_cpue),
        log1p(mean_hh_artisanal_fishing_income),
        use = "complete.obs",
        method = "pearson"
      ),
      spearman_cpue_mean_income = cor(
        mean_cpue,
        mean_hh_artisanal_fishing_income,
        use = "complete.obs",
        method = "spearman"
      ),
      pearson_log_cpue_median_income = cor(
        log1p(mean_cpue),
        log1p(median_hh_artisanal_fishing_income),
        use = "complete.obs",
        method = "pearson"
      ),
      spearman_cpue_median_income = cor(
        mean_cpue,
        median_hh_artisanal_fishing_income,
        use = "complete.obs",
        method = "spearman"
      )
    )

  cpue_income_correlations %>%
    mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
    kable(caption = "Simple correlations between municipality-year CPUE and HHS-derived artisanal-fishing income")
}
Simple correlations between municipality-year CPUE and HHS-derived artisanal-fishing income
n_municipality_years pearson_log_cpue_mean_income spearman_cpue_mean_income pearson_log_cpue_median_income spearman_cpue_median_income
19 0.304 0.125 0.339 0.207

Interpretation note for this section

The CPUE-income analysis should be treated as a municipality-year exploratory diagnostic. A clear positive association would be consistent with the idea that higher fishing productivity is reflected in higher household income from artisanal fishing. A weak or non-significant association would not necessarily mean CPUE is irrelevant to livelihoods. It could reflect limited catch records, imprecise CPUE estimates, mismatch between managed-access areas and HHS sampling areas, price variation, differences in effort, gear and species composition, market access, or the fact that HHS income is household-level while CPUE is trip-level.