Purpose

This short analysis explores household income composition in the Mozambique HHS dataset. The goal is to understand whether household livelihood composition is associated with stronger or weaker outcomes.

The analysis asks two related questions:

  1. Number of income sources: are households with more income sources associated with higher total real monthly household income and better Sustainable Livelihoods outcomes?
  2. Type/combination of income sources: are some income-source combinations associated with higher or lower total real monthly household income and better or worse Sustainable Livelihoods outcomes?

This is an exploratory association analysis, not a causal impact evaluation. Income-source composition is not randomly assigned, and it may reflect deeper differences in geography, market access, wealth, assets, fishery dependence, and program exposure. For that reason, the statistical models adjust for survey year and municipality, and results are interpreted as diagnostic patterns rather than causal effects.

Data and methods

## Read data
# Main HHS file.
# Use the updated dataset with income-source composition if available.
hhs_file <- here("data", "raw", "all_hhs_moz(1).csv")

if (!file.exists(hhs_file)) {
  hhs_file <- here("data", "raw", "all_hhs_moz.csv")
}

hhs_raw <- read_csv(hhs_file, show_col_types = FALSE) %>%
  clean_names()

hhs_raw %>%
  summarise(
    records = n(),
    years = n_distinct(year, na.rm = TRUE),
    municipalities = n_distinct(g1_municipality, na.rm = TRUE),
    communities = n_distinct(g1_community, na.rm = TRUE)
  ) %>%
  kable(caption = "HHS records available for the income-composition analysis")
HHS records available for the income-composition analysis
records years municipalities communities
7297 6 9 27

Income conversion and income-source preparation

The HHS income variable is treated as nominal average monthly household income in MZN. The analysis converts it to constant 2021 MZN using the CPI deflators used in the previous Mozambique income analysis. Values for 2025 and 2026 should be treated as provisional until final annual CPI values are available.

For income-composition analysis, each source is interpreted as the share of household income reported for that source. Missing values in individual source columns are treated as 0% from that source. A source is considered an active income source if it contributes at least 10% of reported household income. This threshold avoids counting very small or possibly incidental sources as meaningful livelihood diversification.

The main analysis excludes records whose reported income-source shares do not sum to roughly 100%. Specifically, it keeps records with total reported source shares between 80% and 120%. This preserves rounding/multiple-response noise while removing clearly inconsistent income-composition records.

# CPI deflators to constant 2021 MZN -------------------------------------

cpi_deflators <- tibble::tribble(
  ~year, ~deflator_to_2021_mzn,
  2019, 1.1012,
  2020, 1.0641,
  2021, 1.0000,
  2022, 0.9068,
  2023, 0.8465,
  2024, 0.8133,
  2025, 0.7792,
  2026, 0.7464
)

# Income-source columns ---------------------------------------------------

source_lookup <- tibble::tribble(
  ~source_col, ~source_label, ~source_short,
  "g4_hh_average_income_source_a_income_farming", "Farming", "farming",
  "g4_hh_average_income_source_b_income_harvesting", "Harvesting", "harvesting",
  "g4_hh_average_income_source_c_income_fishing_artisanal", "Artisanal fishing", "artisanal_fishing",
  "g4_hh_average_income_source_d_income_fishing_aquaculture", "Aquaculture", "aquaculture",
  "g4_hh_average_income_source_e_income_buying_trading", "Buying/trading", "buying_trading",
  "g4_hh_average_income_source_f_income_processing", "Processing", "processing",
  "g4_hh_average_income_source_g_income_extraction", "Extraction", "extraction",
  "g4_hh_average_income_source_h_income_tourism", "Tourism", "tourism",
  "g4_hh_average_income_source_i_income_other_wage", "Other wage", "other_wage",
  "g4_hh_average_income_source_j_income_industrial", "Industrial fishing", "industrial_fishing",
  "g4_hh_average_income_source_k_income_other", "Other", "other"
)

source_cols <- source_lookup$source_col

active_source_threshold <- 10
valid_source_sum_low <- 80
valid_source_sum_high <- 120

# Helper to parse source percentages -------------------------------------

parse_source_pct <- function(x) {
  out <- readr::parse_number(as.character(x))
  out <- dplyr::if_else(is.na(out), 0, out)
  out <- dplyr::if_else(out < 0 | out > 100, NA_real_, out)
  out
}

hhs_income <- hhs_raw %>%
  mutate(
    year = as.integer(year),
    nominal_monthly_hh_income_mzn = readr::parse_number(as.character(g13_hh_average_income))
  ) %>%
  left_join(cpi_deflators, by = "year") %>%
  mutate(
    real_monthly_hh_income_2021_mzn = nominal_monthly_hh_income_mzn * deflator_to_2021_mzn,
    food_security = case_when(
      str_to_lower(g11_food_worry) == "never" ~ 1,
      str_to_lower(g11_food_worry) %in% c("sometimes", "often") ~ 0,
      TRUE ~ NA_real_
    ),
    income_sufficiency = case_when(
      g13_hh_ends_meet %in% c("Fairly easy", "Easy", "Very easy") ~ 1,
      g13_hh_ends_meet %in% c("With difficulty", "With great difficulty") ~ 0,
      TRUE ~ NA_real_
    )
  ) %>%
  mutate(
    across(
      all_of(source_cols),
      parse_source_pct,
      .names = "{.col}_pct"
    )
  )

source_pct_cols <- paste0(source_cols, "_pct")

hhs_income <- hhs_income %>%
  mutate(
    income_source_sum = rowSums(across(all_of(source_pct_cols)), na.rm = TRUE),
    source_composition_valid = between(
      income_source_sum,
      valid_source_sum_low,
      valid_source_sum_high
    )
  ) %>%
  mutate(
    farming_pct = g4_hh_average_income_source_a_income_farming_pct,
    artisanal_fishing_pct = g4_hh_average_income_source_c_income_fishing_artisanal_pct,
    harvesting_pct = g4_hh_average_income_source_b_income_harvesting_pct,
    aquaculture_pct = g4_hh_average_income_source_d_income_fishing_aquaculture_pct,
    buying_trading_pct = g4_hh_average_income_source_e_income_buying_trading_pct,
    processing_pct = g4_hh_average_income_source_f_income_processing_pct,
    extraction_pct = g4_hh_average_income_source_g_income_extraction_pct,
    tourism_pct = g4_hh_average_income_source_h_income_tourism_pct,
    other_wage_pct = g4_hh_average_income_source_i_income_other_wage_pct,
    industrial_fishing_pct = g4_hh_average_income_source_j_income_industrial_pct,
    other_pct = g4_hh_average_income_source_k_income_other_pct
  ) %>%
  mutate(
    active_farming = farming_pct >= active_source_threshold,
    active_harvesting = harvesting_pct >= active_source_threshold,
    active_artisanal_fishing = artisanal_fishing_pct >= active_source_threshold,
    active_aquaculture = aquaculture_pct >= active_source_threshold,
    active_buying_trading = buying_trading_pct >= active_source_threshold,
    active_processing = processing_pct >= active_source_threshold,
    active_extraction = extraction_pct >= active_source_threshold,
    active_tourism = tourism_pct >= active_source_threshold,
    active_other_wage = other_wage_pct >= active_source_threshold,
    active_industrial_fishing = industrial_fishing_pct >= active_source_threshold,
    active_other = other_pct >= active_source_threshold,
    n_active_sources = rowSums(
      across(starts_with("active_")),
      na.rm = TRUE
    ),
    n_active_sources_cat = case_when(
      n_active_sources >= 4 ~ "4+",
      TRUE ~ as.character(n_active_sources)
    ),
    n_active_sources_cat = factor(
      n_active_sources_cat,
      levels = c("0", "1", "2", "3", "4+")
    )
  ) %>%
  rowwise() %>%
  mutate(
    income_combination = {
      active_sources <- c(
        "Farming" = active_farming,
        "Harvesting" = active_harvesting,
        "Artisanal fishing" = active_artisanal_fishing,
        "Aquaculture" = active_aquaculture,
        "Buying/trading" = active_buying_trading,
        "Processing" = active_processing,
        "Extraction" = active_extraction,
        "Tourism" = active_tourism,
        "Other wage" = active_other_wage,
        "Industrial fishing" = active_industrial_fishing,
        "Other" = active_other
      )
      out <- names(active_sources)[active_sources]
      if (length(out) == 0) "No source >=10%" else paste(out, collapse = " + ")
    }
  ) %>%
  ungroup()

# Clean income for analysis ----------------------------------------------
# For total-income analysis we remove extreme values above the year-specific p99.
# This is applied only to the total-income models and plots.

hhs_income <- hhs_income %>%
  group_by(year) %>%
  mutate(
    income_p99_year = quantile(real_monthly_hh_income_2021_mzn, 0.99, na.rm = TRUE)
  ) %>%
  ungroup() %>%
  mutate(
    income_clean = !is.na(real_monthly_hh_income_2021_mzn) &
      real_monthly_hh_income_2021_mzn >= 0 &
      real_monthly_hh_income_2021_mzn <= income_p99_year,
    log1p_real_income = log1p(real_monthly_hh_income_2021_mzn)
  )

composition_data <- hhs_income %>%
  filter(source_composition_valid)

model_income_data <- composition_data %>%
  filter(income_clean)

composition_coverage <- hhs_income %>%
  summarise(
    total_records = n(),
    valid_source_composition = sum(source_composition_valid, na.rm = TRUE),
    excluded_source_composition = sum(!source_composition_valid, na.rm = TRUE),
    income_records_after_p99_cleaning = sum(income_clean & source_composition_valid, na.rm = TRUE)
  )

composition_coverage %>%
  kable(caption = "Records retained for income-composition analysis")
Records retained for income-composition analysis
total_records valid_source_composition excluded_source_composition income_records_after_p99_cleaning
7297 6632 665 6277

Modeling approach

The statistical analysis uses three complementary models:

  1. Total real household income: linear models of log(1 + real monthly household income in 2021 MZN). Effects are presented as approximate percent differences in total real monthly household income.
  2. Income sufficiency / ability to meet needs: logistic models where 1 means the household reports covering needs fairly easily, easily, or very easily.
  3. Food security: logistic models where 1 means the household never worried about not having enough food in the last 12 months.

All models adjust for survey year and municipality. For source-combination models, combinations with very small samples are excluded from the model ranking to avoid unstable estimates. P-values for multiple source-combination comparisons are adjusted using the Benjamini-Hochberg false-discovery-rate procedure.

# Summary helpers ---------------------------------------------------------

prop_ci <- function(x) {
  x <- x[!is.na(x)]
  n <- length(x)
  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)
  )
}

tidy_wald <- function(model) {
  broom::tidy(model) %>%
    mutate(
      conf.low = estimate - 1.96 * std.error,
      conf.high = estimate + 1.96 * std.error
    )
}

extract_lm_factor_pct <- function(model, factor_prefix) {
  tidy_wald(model) %>%
    filter(str_starts(term, factor_prefix)) %>%
    mutate(
      level = str_remove(term, paste0("^", factor_prefix)),
      percent_difference = 100 * (exp(estimate) - 1),
      ci_low = 100 * (exp(conf.low) - 1),
      ci_high = 100 * (exp(conf.high) - 1),
      p_value = p.value
    ) %>%
    select(level, percent_difference, ci_low, ci_high, p_value)
}

extract_glm_factor_or <- function(model, factor_prefix) {
  tidy_wald(model) %>%
    filter(str_starts(term, factor_prefix)) %>%
    mutate(
      level = str_remove(term, paste0("^", factor_prefix)),
      odds_ratio = exp(estimate),
      ci_low = exp(conf.low),
      ci_high = exp(conf.high),
      p_value = p.value
    ) %>%
    select(level, odds_ratio, ci_low, ci_high, p_value)
}

fmt_p <- function(p) {
  case_when(
    is.na(p) ~ NA_character_,
    p < 0.001 ~ "<0.001",
    TRUE ~ sprintf("%.3f", p)
  )
}

1. Income composition overview

Average income-source shares by year

This figure shows how the average reported income composition changes through time across the Mozambique HHS sample. Because the HHS coverage changes across municipalities and years, this plot should be interpreted as a descriptive overview of the available sample rather than a nationally representative trend.

source_long <- composition_data %>%
  select(year, g1_municipality, all_of(source_pct_cols)) %>%
  pivot_longer(
    cols = all_of(source_pct_cols),
    names_to = "source_col_pct",
    values_to = "source_pct"
  ) %>%
  mutate(
    source_col = str_remove(source_col_pct, "_pct$")
  ) %>%
  left_join(source_lookup, by = "source_col")

source_year <- source_long %>%
  group_by(year, source_label) %>%
  summarise(
    mean_share = mean(source_pct, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(
  source_year,
  aes(
    x = year,
    y = mean_share,
    fill = source_label
  )
) +
  geom_area(alpha = 0.9, color = "white", linewidth = 0.2) +
  scale_y_continuous(labels = label_percent(scale = 1)) +
  scale_x_continuous(breaks = sort(unique(source_year$year))) +
  labs(
    title = "Average reported household income composition by year",
    subtitle = "Source shares are based on records where income-source percentages sum to roughly 100%",
    x = "Survey year",
    y = "Average share of household income",
    fill = "Income source"
  ) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Municipality-level source composition

source_municipality <- source_long %>%
  group_by(g1_municipality, source_label) %>%
  summarise(
    mean_share = mean(source_pct, na.rm = TRUE),
    .groups = "drop"
  )

source_municipality %>%
  mutate(
    g1_municipality = fct_reorder(g1_municipality, mean_share, .fun = sum),
    label = if_else(mean_share >= 5, paste0(round(mean_share, 0), "%"), ""),
    label_color = if_else(mean_share > 25, "white", "black")
  ) %>%
  ggplot(
    aes(
      x = source_label,
      y = g1_municipality,
      fill = mean_share
    )
  ) +
  geom_tile(color = "white") +
  geom_text(aes(label = label, color = label_color), size = 3) +
  scale_color_identity() +
  scale_fill_viridis_c(labels = label_percent(scale = 1), option = "C") +
  labs(
    title = "Average income-source shares by municipality",
    subtitle = "Cell values shown where the average source share is at least 5%",
    x = NULL,
    y = NULL,
    fill = "Average share"
  ) +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

2. Number of income sources

Distribution of number of active income sources

An active income source is defined as a source contributing at least 10% of household income.

composition_data %>%
  count(n_active_sources_cat) %>%
  mutate(
    pct = n / sum(n)
  ) %>%
  ggplot(aes(x = n_active_sources_cat, y = pct)) +
  geom_col(width = 0.7) +
  geom_text(aes(label = percent(pct, accuracy = 0.1)), vjust = -0.4) +
  scale_y_continuous(labels = label_percent(), limits = c(0, NA)) +
  labs(
    title = "Number of active household income sources",
    subtitle = paste0("Active source = at least ", active_source_threshold, "% of household income"),
    x = "Number of active income sources",
    y = "Share of HHS records"
  )

composition_data %>%
  count(g1_municipality, n_active_sources_cat) %>%
  group_by(g1_municipality) %>%
  mutate(pct = n / sum(n)) %>%
  ungroup() %>%
  ggplot(aes(x = n_active_sources_cat, y = pct, fill = n_active_sources_cat)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ g1_municipality, ncol = 3) +
  scale_y_continuous(labels = label_percent()) +
  labs(
    title = "Number of active income sources by municipality",
    x = "Number of active sources",
    y = "Share of HHS records"
  ) +
  theme(panel.grid.minor = element_blank())

Number of sources and total real household income

income_by_n_sources <- model_income_data %>%
  group_by(n_active_sources_cat) %>%
  summarise(
    n = n(),
    median_income = median(real_monthly_hh_income_2021_mzn, na.rm = TRUE),
    q25 = quantile(real_monthly_hh_income_2021_mzn, 0.25, na.rm = TRUE),
    q75 = quantile(real_monthly_hh_income_2021_mzn, 0.75, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(
  income_by_n_sources,
  aes(x = n_active_sources_cat, y = median_income)
) +
  geom_col(width = 0.7, fill = "grey55") +
  geom_errorbar(aes(ymin = q25, ymax = q75), width = 0.2) +
  geom_text(aes(label = paste0("n=", n)), vjust = -0.4, size = 3) +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Median real household income by number of active income sources",
    subtitle = "Bars show median; error bars show IQR; income is constant 2021 MZN",
    x = "Number of active income sources",
    y = "Median monthly household income, 2021 MZN"
  )

model_income_data %>%
  group_by(g1_municipality, n_active_sources_cat) %>%
  summarise(
    n = n(),
    median_income = median(real_monthly_hh_income_2021_mzn, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = n_active_sources_cat, y = median_income)) +
  geom_col(width = 0.7, fill = "grey55") +
  geom_text(aes(label = paste0("n=", n)), vjust = -0.4, size = 2.8) +
  facet_wrap(~ g1_municipality, ncol = 3, scales = "free_y") +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Median real household income by number of active sources and municipality",
    x = "Number of active income sources",
    y = "Median monthly household income, 2021 MZN"
  ) +
  theme(panel.grid.minor = element_blank())

Adjusted model: number of sources and total income

The model below estimates the association between the number of active income sources and total real household income, adjusting for survey year and municipality. The reference category is households with one active income source.

income_n_sources_model_data <- model_income_data %>%
  mutate(
    n_active_sources_ref = relevel(n_active_sources_cat, ref = "1")
  )

model_income_n_sources <- lm(
  log1p_real_income ~ n_active_sources_ref + factor(year) + factor(g1_municipality),
  data = income_n_sources_model_data
)

income_n_sources_results <- extract_lm_factor_pct(
  model_income_n_sources,
  "n_active_sources_ref"
) %>%
  mutate(
    level = recode(level, "4+" = "4+"),
    p_value_label = fmt_p(p_value)
  )

income_n_sources_results %>%
  mutate(
    across(c(percent_difference, ci_low, ci_high), ~ round(.x, 1))
  ) %>%
  kable(
    caption = "Adjusted association between number of active income sources and total real household income"
  )
Adjusted association between number of active income sources and total real household income
level percent_difference ci_low ci_high p_value p_value_label
2 28.4 18.0 39.7 0.0000000 <0.001
3 34.5 18.0 53.4 0.0000096 <0.001
4+ 30.5 11.9 52.2 0.0006792 <0.001 
income_n_sources_results %>%
  ggplot(aes(x = percent_difference, y = fct_reorder(level, percent_difference))) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey40") +
  geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.2, color = "grey35") +
  geom_point(size = 2.6) +
  labs(
    title = "Adjusted association: number of active income sources and total income",
    subtitle = "Percent difference in total real monthly household income relative to one source; adjusted for year and municipality",
    x = "% difference in total real monthly household income",
    y = "Number of active income sources"
  ) +
  theme(panel.grid.minor = element_blank())

The main takeaway is that households with more than one active income source have significantly higher total real monthly household income than households with only one active income source, even after adjusting for year and municipality.

More specifically, compared with households with one source, households with:

  • 2 sources have about 28% higher real monthly income.

  • 3 sources have about 35% higher real monthly income.

  • 4+ sources have about 31% higher real monthly income.

All of these associations are statistically significant, and the confidence intervals are clearly above zero. The pattern supports the livelihood-diversification idea: having multiple meaningful income sources is associated with higher household income.

3. Income-source combinations and total income

Common income-source combinations

The source-combination analysis focuses on combinations with enough records to support interpretation. By default, this includes combinations with at least 50 HHS records after income cleaning and valid source-composition filtering.

min_combo_n <- 50

combo_counts <- model_income_data %>%
  count(income_combination, n_active_sources, sort = TRUE) %>%
  mutate(
    include_in_combo_model = n >= min_combo_n
  )

combo_counts %>%
  slice_head(n = 25) %>%
  kable(caption = "Most common income-source combinations")
Most common income-source combinations
income_combination n_active_sources n include_in_combo_model
Farming + Artisanal fishing 2 1356 TRUE
Artisanal fishing 1 961 TRUE
Buying/trading 1 582 TRUE
Farming + Buying/trading 2 500 TRUE
Farming 1 340 TRUE
Farming + Artisanal fishing + Buying/trading 3 330 TRUE
Farming + Other 2 182 TRUE
Other 1 177 TRUE
Farming + Artisanal fishing + Buying/trading + Processing 4 169 TRUE
Buying/trading + Processing 2 128 TRUE
Artisanal fishing + Aquaculture + Buying/trading + Industrial fishing 4 112 TRUE
Artisanal fishing + Other 2 96 TRUE
Artisanal fishing + Buying/trading 2 87 TRUE
Farming + Buying/trading + Processing 3 78 TRUE
Farming + Artisanal fishing + Other 3 77 TRUE
Farming + Artisanal fishing + Processing 3 73 TRUE
Processing 1 64 TRUE
Farming + Harvesting 2 58 TRUE
Other wage 1 51 TRUE
Artisanal fishing + Industrial fishing 2 49 FALSE
Artisanal fishing + Processing 2 45 FALSE
Farming + Harvesting + Artisanal fishing + Buying/trading + Processing 5 44 FALSE
Farming + Other wage 2 39 FALSE
Artisanal fishing + Buying/trading + Processing 3 37 FALSE
Farming + Harvesting + Artisanal fishing 3 36 FALSE 
combo_counts %>%
  slice_head(n = 20) %>%
  mutate(
    income_combination = fct_reorder(str_wrap(income_combination, 45), n)
  ) %>%
  ggplot(aes(x = n, y = income_combination)) +
  geom_col(fill = "grey55") +
  labs(
    title = "Most common income-source combinations",
    subtitle = paste0("Active source = at least ", active_source_threshold, "% of household income"),
    x = "HHS records",
    y = NULL
  )

Median income by common income-source combination

combo_income_summary <- model_income_data %>%
  semi_join(
    combo_counts %>% filter(n >= min_combo_n),
    by = c("income_combination", "n_active_sources")
  ) %>%
  group_by(income_combination, n_active_sources) %>%
  summarise(
    n = n(),
    median_income = median(real_monthly_hh_income_2021_mzn, na.rm = TRUE),
    q25 = quantile(real_monthly_hh_income_2021_mzn, 0.25, na.rm = TRUE),
    q75 = quantile(real_monthly_hh_income_2021_mzn, 0.75, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(
  combo_income_summary %>%
    mutate(income_combination = fct_reorder(str_wrap(income_combination, 45), median_income)),
  aes(x = median_income, y = income_combination)
) +
  geom_errorbarh(aes(xmin = q25, xmax = q75), height = 0.2, color = "grey60") +
  geom_point(aes(size = n), alpha = 0.8) +
  scale_x_continuous(labels = comma) +
  labs(
    title = "Median real household income by income-source combination",
    subtitle = "Points show median; horizontal bars show IQR; point size shows number of HHS records",
    x = "Median monthly household income, 2021 MZN",
    y = NULL,
    size = "N"
  ) +
  theme(panel.grid.minor = element_blank())

Adjusted model: source combinations and total income

The model below compares common source combinations to the most common combination in the data. This gives a statistically justified answer to whether each common combination is associated with higher or lower total income than the reference combination, after adjusting for year and municipality. P-values are FDR-adjusted because many combinations are compared at the same time.

combo_reference <- combo_counts %>%
  filter(n >= min_combo_n) %>%
  arrange(desc(n)) %>%
  slice(1) %>%
  pull(income_combination)

combo_model_data <- model_income_data %>%
  semi_join(
    combo_counts %>% filter(n >= min_combo_n),
    by = c("income_combination", "n_active_sources")
  ) %>%
  mutate(
    combo_factor = factor(income_combination),
    combo_factor = relevel(combo_factor, ref = combo_reference)
  )

model_income_combo <- lm(
  log1p_real_income ~ combo_factor + factor(year) + factor(g1_municipality),
  data = combo_model_data
)

income_combo_results <- extract_lm_factor_pct(model_income_combo, "combo_factor") %>%
  mutate(
    fdr_p_value = p.adjust(p_value, method = "BH"),
    p_value_label = fmt_p(p_value),
    fdr_p_value_label = fmt_p(fdr_p_value),
    direction = case_when(
      fdr_p_value < 0.05 & percent_difference > 0 ~ "Statistically higher than reference",
      fdr_p_value < 0.05 & percent_difference < 0 ~ "Statistically lower than reference",
      TRUE ~ "Not statistically clear"
    )
  )

combo_reference
## [1] "Farming + Artisanal fishing"
income_combo_results %>%
  arrange(percent_difference) %>%
  mutate(
    across(c(percent_difference, ci_low, ci_high), ~ round(.x, 1))
  ) %>%
  select(level, percent_difference, ci_low, ci_high, p_value_label, fdr_p_value_label, direction) %>%
  kable(
    caption = paste0("Adjusted total-income differences by income-source combination; reference = ", combo_reference)
  )
Adjusted total-income differences by income-source combination; reference = Farming + Artisanal fishing
level percent_difference ci_low ci_high p_value_label fdr_p_value_label direction
Farming -74.7 -78.6 -70.1 <0.001 <0.001 Statistically lower than reference
Farming + Harvesting -52.5 -67.0 -31.5 <0.001 <0.001 Statistically lower than reference
Farming + Buying/trading -11.4 -23.5 2.6 0.107 0.148 Not statistically clear
Farming + Other -11.1 -28.2 10.0 0.278 0.334 Not statistically clear
Farming + Artisanal fishing + Buying/trading + Processing -6.2 -27.2 21.0 0.623 0.693 Not statistically clear
Other -5.1 -24.6 19.4 0.654 0.693 Not statistically clear
Farming + Artisanal fishing + Processing 2.2 -26.5 42.1 0.896 0.896 Not statistically clear
Artisanal fishing 8.4 -3.9 22.4 0.189 0.244 Not statistically clear
Farming + Artisanal fishing + Buying/trading 16.4 -2.7 39.2 0.098 0.147 Not statistically clear
Buying/trading 21.5 5.6 39.9 0.006 0.019 Statistically higher than reference
Processing 35.3 -3.7 90.2 0.082 0.134 Not statistically clear
Artisanal fishing + Buying/trading 37.1 2.1 84.1 0.036 0.072 Not statistically clear
Farming + Artisanal fishing + Other 39.4 2.2 90.2 0.036 0.072 Not statistically clear
Other wage 48.9 1.3 119.0 0.043 0.077 Not statistically clear
Farming + Buying/trading + Processing 51.1 8.1 111.2 0.016 0.040 Statistically higher than reference
Artisanal fishing + Other 54.8 16.7 105.2 0.002 0.009 Statistically higher than reference
Artisanal fishing + Aquaculture + Buying/trading + Industrial fishing 120.1 67.4 189.3 <0.001 <0.001 Statistically higher than reference
Buying/trading + Processing 167.6 107.6 245.0 <0.001 <0.001 Statistically higher than reference 
income_combo_results %>%
  mutate(
    level = fct_reorder(str_wrap(level, 45), percent_difference)
  ) %>%
  ggplot(aes(x = percent_difference, y = level, color = direction)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey40") +
  geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.2, alpha = 0.7) +
  geom_point(size = 2.6) +
  labs(
    title = "Adjusted association: source combinations and total real household income",
    subtitle = paste0("Reference combination = ", combo_reference, "; adjusted for year and municipality; FDR used for multiple comparisons"),
    x = "% difference in total real monthly household income",
    y = NULL,
    color = NULL
  ) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Using Farming + Artisanal fishing as the reference group, the combinations that look statistically higher are mainly those involving buying/trading and/or processing, especially Buying/trading + Processing, which has the largest positive association. Some mixed combinations that include artisanal fishing, buying/trading, aquaculture, or industrial fishing also appear higher, though some have wider uncertainty.

In contrast, Farming only is strongly and statistically associated with lower total household income, and Farming + Harvesting is also lower. Many other combinations are not statistically distinguishable from the reference group after adjustment.

4. Number of income sources and Sustainable Livelihoods

This section repeats the analysis using the two Sustainable Livelihoods indicators:

  • Food security: respondent never worried about having enough food for everyone in the household.
  • Income sufficiency / ability to meet needs: household can cover its needs fairly easily, easily, or very easily.

Descriptive patterns

sl_by_n_sources <- composition_data %>%
  select(n_active_sources_cat, food_security, income_sufficiency) %>%
  pivot_longer(
    cols = c(food_security, income_sufficiency),
    names_to = "outcome",
    values_to = "value"
  ) %>%
  mutate(
    outcome = recode(
      outcome,
      food_security = "Food security",
      income_sufficiency = "Income sufficiency / ability to meet needs"
    )
  ) %>%
  group_by(n_active_sources_cat, outcome) %>%
  summarise(prop_ci(value), .groups = "drop")

ggplot(
  sl_by_n_sources,
  aes(
    x = n_active_sources_cat,
    y = pct,
    color = outcome,
    group = outcome
  )
) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.15, alpha = 0.7) +
  geom_line(linewidth = 0.8) +
  geom_point(size = 2.5) +
  scale_y_continuous(labels = label_percent(scale = 1), limits = c(0, 100)) +
  labs(
    title = "Sustainable Livelihoods indicators by number of active income sources",
    subtitle = "Error bars show approximate 95% CIs for the proportion",
    x = "Number of active income sources",
    y = "% positive",
    color = "Indicator"
  ) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Adjusted models: number of sources and Sustainable Livelihoods

The logistic models below estimate whether households with different numbers of active income sources have higher or lower odds of positive Sustainable Livelihoods outcomes, relative to households with one active income source, adjusting for year and municipality.

sl_n_sources_model_data <- composition_data %>%
  mutate(
    n_active_sources_ref = relevel(n_active_sources_cat, ref = "1")
  )

model_food_n_sources <- glm(
  food_security ~ n_active_sources_ref + factor(year) + factor(g1_municipality),
  data = sl_n_sources_model_data,
  family = binomial()
)

model_needs_n_sources <- glm(
  income_sufficiency ~ n_active_sources_ref + factor(year) + factor(g1_municipality),
  data = sl_n_sources_model_data,
  family = binomial()
)

sl_n_sources_results <- bind_rows(
  extract_glm_factor_or(model_food_n_sources, "n_active_sources_ref") %>%
    mutate(outcome = "Food security"),
  extract_glm_factor_or(model_needs_n_sources, "n_active_sources_ref") %>%
    mutate(outcome = "Income sufficiency / ability to meet needs")
) %>%
  mutate(
    p_value_label = fmt_p(p_value)
  )

sl_n_sources_results %>%
  mutate(
    across(c(odds_ratio, ci_low, ci_high), ~ round(.x, 3))
  ) %>%
  select(outcome, level, odds_ratio, ci_low, ci_high, p_value_label) %>%
  kable(caption = "Adjusted association between number of active sources and Sustainable Livelihoods outcomes")
Adjusted association between number of active sources and Sustainable Livelihoods outcomes
outcome level odds_ratio ci_low ci_high p_value_label
Food security 2 0.962 0.792 1.169 0.699
Food security 3 0.486 0.349 0.677 <0.001
Food security 4+ 0.118 0.073 0.193 <0.001
Income sufficiency / ability to meet needs 2 0.659 0.575 0.756 <0.001
Income sufficiency / ability to meet needs 3 0.517 0.412 0.649 <0.001
Income sufficiency / ability to meet needs 4+ 0.448 0.345 0.584 <0.001 
sl_n_sources_results %>%
  mutate(
    level = fct_reorder(level, odds_ratio)
  ) %>%
  ggplot(aes(x = odds_ratio, y = level)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey40") +
  geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.2, color = "grey35") +
  geom_point(size = 2.6) +
  facet_wrap(~ outcome, ncol = 1) +
  scale_x_log10() +
  labs(
    title = "Adjusted association: number of sources and Sustainable Livelihoods outcomes",
    subtitle = "Odds ratios relative to households with one active income source; adjusted for year and municipality",
    x = "Odds ratio, log scale",
    y = "Number of active income sources"
  ) +
  theme(panel.grid.minor = element_blank())

Having more active income sources is associated with higher total income, but not necessarily with better Sustainable Livelihoods outcomes.

In fact, after adjusting for year and municipality, households with more income sources have lower odds of reporting food security and income sufficiency compared with households with only one active income source. The pattern is especially strong for 3 or 4+ sources, where the odds of positive livelihood outcomes are clearly below 1.

5. Source combinations and Sustainable Livelihoods

Descriptive Sustainable Livelihoods by common combinations

sl_combo_summary <- composition_data %>%
  semi_join(
    combo_counts %>% filter(n >= min_combo_n),
    by = c("income_combination", "n_active_sources")
  ) %>%
  select(income_combination, food_security, income_sufficiency) %>%
  pivot_longer(
    cols = c(food_security, income_sufficiency),
    names_to = "outcome",
    values_to = "value"
  ) %>%
  mutate(
    outcome = recode(
      outcome,
      food_security = "Food security",
      income_sufficiency = "Income sufficiency / ability to meet needs"
    )
  ) %>%
  group_by(income_combination, outcome) %>%
  summarise(prop_ci(value), .groups = "drop")

ggplot(
  sl_combo_summary %>%
    mutate(income_combination = fct_reorder(str_wrap(income_combination, 45), pct)),
  aes(x = pct, y = income_combination, color = outcome)
) +
  geom_errorbarh(aes(xmin = ci_low, xmax = ci_high), height = 0.2, alpha = 0.6) +
  geom_point(size = 2.2) +
  scale_x_continuous(labels = label_percent(scale = 1), limits = c(0, 100)) +
  labs(
    title = "Sustainable Livelihoods indicators by common income-source combination",
    subtitle = "Points show proportion positive; horizontal bars show approximate 95% CIs",
    x = "% positive",
    y = NULL,
    color = "Indicator"
  ) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Adjusted models: combinations and Sustainable Livelihoods

The models below compare common income-source combinations to the most common reference combination, adjusting for year and municipality. P-values are FDR-adjusted across the source-combination comparisons within each outcome.

combo_sl_model_data <- composition_data %>%
  semi_join(
    combo_counts %>% filter(n >= min_combo_n),
    by = c("income_combination", "n_active_sources")
  ) %>%
  mutate(
    combo_factor = factor(income_combination),
    combo_factor = relevel(combo_factor, ref = combo_reference)
  )

model_food_combo <- glm(
  food_security ~ combo_factor + factor(year) + factor(g1_municipality),
  data = combo_sl_model_data,
  family = binomial()
)

model_needs_combo <- glm(
  income_sufficiency ~ combo_factor + factor(year) + factor(g1_municipality),
  data = combo_sl_model_data,
  family = binomial()
)

sl_combo_results <- bind_rows(
  extract_glm_factor_or(model_food_combo, "combo_factor") %>%
    mutate(outcome = "Food security"),
  extract_glm_factor_or(model_needs_combo, "combo_factor") %>%
    mutate(outcome = "Income sufficiency / ability to meet needs")
) %>%
  group_by(outcome) %>%
  mutate(
    fdr_p_value = p.adjust(p_value, method = "BH")
  ) %>%
  ungroup() %>%
  mutate(
    p_value_label = fmt_p(p_value),
    fdr_p_value_label = fmt_p(fdr_p_value),
    direction = case_when(
      fdr_p_value < 0.05 & odds_ratio > 1 ~ "Statistically higher odds than reference",
      fdr_p_value < 0.05 & odds_ratio < 1 ~ "Statistically lower odds than reference",
      TRUE ~ "Not statistically clear"
    )
  )

sl_combo_results %>%
  mutate(
    across(c(odds_ratio, ci_low, ci_high), ~ round(.x, 3))
  ) %>%
  select(outcome, level, odds_ratio, ci_low, ci_high, p_value_label, fdr_p_value_label, direction) %>%
  arrange(outcome, odds_ratio) %>%
  kable(
    caption = paste0("Adjusted Sustainable Livelihoods associations by source combination; reference = ", combo_reference)
  )
Adjusted Sustainable Livelihoods associations by source combination; reference = Farming + Artisanal fishing
outcome level odds_ratio ci_low ci_high p_value_label fdr_p_value_label direction
Food security Artisanal fishing + Aquaculture + Buying/trading + Industrial fishing 0.000 0.000 Inf 0.971 0.972 Not statistically clear
Food security Farming + Artisanal fishing + Processing 0.000 0.000 Inf 0.969 0.972 Not statistically clear
Food security Farming + Artisanal fishing + Buying/trading 0.053 0.023 1.230000e-01 <0.001 <0.001 Statistically lower odds than reference
Food security Farming + Artisanal fishing + Buying/trading + Processing 0.112 0.053 2.340000e-01 <0.001 <0.001 Statistically lower odds than reference
Food security Processing 0.378 0.070 2.048000e+00 0.259 0.389 Not statistically clear
Food security Farming 0.558 0.357 8.710000e-01 0.010 0.026 Statistically lower odds than reference
Food security Artisanal fishing + Buying/trading 0.769 0.298 1.983000e+00 0.587 0.813 Not statistically clear
Food security Farming + Harvesting 0.982 0.361 2.670000e+00 0.972 0.972 Not statistically clear
Food security Artisanal fishing 1.022 0.758 1.376000e+00 0.888 0.972 Not statistically clear
Food security Buying/trading + Processing 1.024 0.481 2.182000e+00 0.950 0.972 Not statistically clear
Food security Farming + Buying/trading 1.238 0.875 1.751000e+00 0.228 0.373 Not statistically clear
Food security Farming + Other 1.773 1.027 3.061000e+00 0.040 0.080 Not statistically clear
Food security Artisanal fishing + Other 1.934 0.988 3.785000e+00 0.054 0.097 Not statistically clear
Food security Buying/trading 2.024 1.455 2.816000e+00 <0.001 <0.001 Statistically higher odds than reference
Food security Farming + Buying/trading + Processing 2.048 1.110 3.779000e+00 0.022 0.049 Statistically higher odds than reference
Food security Farming + Artisanal fishing + Other 2.505 1.270 4.941000e+00 0.008 0.024 Statistically higher odds than reference
Food security Other 2.549 1.558 4.170000e+00 <0.001 <0.001 Statistically higher odds than reference
Food security Other wage 8.045 4.151 1.559300e+01 <0.001 <0.001 Statistically higher odds than reference
Income sufficiency / ability to meet needs Artisanal fishing + Aquaculture + Buying/trading + Industrial fishing 0.000 0.000 7.018051e+184 0.947 0.947 Not statistically clear
Income sufficiency / ability to meet needs Farming + Artisanal fishing + Buying/trading 0.257 0.164 4.030000e-01 <0.001 <0.001 Statistically lower odds than reference
Income sufficiency / ability to meet needs Farming + Artisanal fishing + Buying/trading + Processing 0.413 0.255 6.680000e-01 <0.001 <0.001 Statistically lower odds than reference
Income sufficiency / ability to meet needs Farming 0.454 0.327 6.290000e-01 <0.001 <0.001 Statistically lower odds than reference
Income sufficiency / ability to meet needs Farming + Harvesting 0.542 0.261 1.124000e+00 0.100 0.149 Not statistically clear
Income sufficiency / ability to meet needs Farming + Artisanal fishing + Processing 0.808 0.430 1.518000e+00 0.508 0.624 Not statistically clear
Income sufficiency / ability to meet needs Farming + Other 0.938 0.631 1.396000e+00 0.753 0.848 Not statistically clear
Income sufficiency / ability to meet needs Farming + Artisanal fishing + Other 1.039 0.602 1.795000e+00 0.890 0.942 Not statistically clear
Income sufficiency / ability to meet needs Artisanal fishing + Buying/trading 1.190 0.700 2.023000e+00 0.520 0.624 Not statistically clear
Income sufficiency / ability to meet needs Farming + Buying/trading 1.394 1.073 1.812000e+00 0.013 0.026 Statistically higher odds than reference
Income sufficiency / ability to meet needs Farming + Buying/trading + Processing 1.454 0.853 2.480000e+00 0.169 0.234 Not statistically clear
Income sufficiency / ability to meet needs Artisanal fishing 1.474 1.196 1.816000e+00 <0.001 <0.001 Statistically higher odds than reference
Income sufficiency / ability to meet needs Artisanal fishing + Other 1.514 0.952 2.408000e+00 0.079 0.130 Not statistically clear
Income sufficiency / ability to meet needs Buying/trading + Processing 1.644 1.064 2.541000e+00 0.025 0.045 Statistically higher odds than reference
Income sufficiency / ability to meet needs Other 2.510 1.767 3.565000e+00 <0.001 <0.001 Statistically higher odds than reference
Income sufficiency / ability to meet needs Buying/trading 2.854 2.258 3.607000e+00 <0.001 <0.001 Statistically higher odds than reference
Income sufficiency / ability to meet needs Processing 4.029 2.391 6.790000e+00 <0.001 <0.001 Statistically higher odds than reference
Income sufficiency / ability to meet needs Other wage 4.501 2.550 7.943000e+00 <0.001 <0.001 Statistically higher odds than reference 
or_plot_min <- 0.01
or_plot_max <- 10

sl_combo_results_plot <- sl_combo_results %>%
  mutate(
    level_wrapped = str_wrap(level, 45),
    level_wrapped = fct_reorder(level_wrapped, odds_ratio),
    
    # Keep actual estimates in the table, but cap displayed CIs in the plot
    ci_low_plot = pmax(ci_low, or_plot_min),
    ci_high_plot = pmin(ci_high, or_plot_max),
    
    # Also cap displayed point if an estimate itself is outside the plot range
    odds_ratio_plot = pmin(pmax(odds_ratio, or_plot_min), or_plot_max),
    
    ci_truncated = ci_low < or_plot_min | ci_high > or_plot_max
  )

ggplot(
  sl_combo_results_plot,
  aes(
    x = odds_ratio_plot,
    y = level_wrapped,
    color = direction
  )
) +
  geom_vline(
    xintercept = 1,
    linetype = "dashed",
    color = "grey40"
  ) +
  geom_segment(
    aes(
      x = ci_low_plot,
      xend = ci_high_plot,
      y = level_wrapped,
      yend = level_wrapped
    ),
    linewidth = 0.6,
    alpha = 0.75
  ) +
  geom_point(size = 2.4) +
  facet_wrap(
    ~ outcome,
    ncol = 1,
    scales = "free_y"
  ) +
  scale_x_log10(
    limits = c(or_plot_min, or_plot_max),
    breaks = c(0.1, 0.25, 0.5, 1, 2, 5, 10),
    labels = c("0.1", "0.25", "0.5", "1", "2", "5", "10")
  ) +
  labs(
    title = "Adjusted association: source combinations and Sustainable Livelihoods outcomes",
    subtitle = paste0(
      "Reference combination = ",
      combo_reference,
      "; adjusted for year and municipality; FDR used for multiple comparisons; displayed CIs truncated at 0.1–10"
    ),
    x = "Odds ratio, log scale",
    y = NULL,
    color = NULL
  ) +
  theme(
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

Interpretation notes

The core interpretation should separate descriptive patterns from adjusted statistical associations. The descriptive plots show which income-source structures appear more common and how outcomes vary across them. The adjusted models provide the statistical basis for saying whether a difference remains after accounting for survey year and municipality.

Suggested cautious framing:

Households with more income sources and some source combinations may show higher total real monthly income or better Sustainable Livelihoods outcomes, but these patterns should be interpreted as associations rather than causal effects. Income-source composition likely reflects underlying differences in opportunity, geography, market access, household assets, and program exposure. The models adjust for year and municipality, but they do not fully control for all household-level differences that shape both livelihood composition and outcomes.

For program interpretation, the most useful finding is not simply whether diversification is always good or bad, but whether particular livelihood portfolios appear consistently associated with stronger or weaker outcomes. Where a source combination is statistically associated with lower income or lower Sustainable Livelihoods outcomes, it may indicate households that are more vulnerable or livelihood strategies that need stronger support, market linkages, or risk reduction.

{r export-tables} # Export key tables # output_dir <- here("outputs") # if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE) # # write_csv(composition_coverage, here("outputs", "income_composition_coverage.csv")) # write_csv(income_n_sources_results, here("outputs", "model_number_sources_total_income.csv")) # write_csv(income_combo_results, here("outputs", "model_combinations_total_income.csv")) # write_csv(sl_n_sources_results, here("outputs", "model_number_sources_sustainable_livelihoods.csv")) # write_csv(sl_combo_results, here("outputs", "model_combinations_sustainable_livelihoods.csv")) # write_csv(combo_counts, here("outputs", "income_combination_counts.csv")) #