# =============================================================================
# KARANGURA CFCT BASELINE
# PAPER 3 ANALYSIS SCRIPT
#
# Manuscript:
# Community Water Management, Environmental Stewardship and Disaster Resilience
# in Karangura, Uganda:
# Evidence from a Rural Household Baseline Survey
#
# Input folder:
#   C:/Users/Stephen/Documents/Researches/Karangura Papers/
#     03_Paper_Water_Governance_Disaster_Resilience
#
# Required input in that folder:
#   Karangura_BMC_master_cleaned_household_dataset.csv
#
# Outputs:
#   Paper3_Outputs_YYYYMMDD/Tables
#   Paper3_Outputs_YYYYMMDD/Figures
#   Paper3_Outputs_YYYYMMDD/Manuscript
#   Paper3_Outputs_YYYYMMDD/R_objects
#   Paper3_Outputs_YYYYMMDD/Logs
# =============================================================================

rm(list = ls())
gc()

required_packages <- c(
  "tidyverse", "janitor", "openxlsx", "flextable", "officer",
  "broom", "sandwich", "lmtest", "scales", "stringr", "forcats"
)
new_packages <- required_packages[!required_packages %in% installed.packages()[, "Package"]]
if (length(new_packages) > 0) install.packages(new_packages, dependencies = TRUE)

suppressPackageStartupMessages({
  library(tidyverse)
  library(janitor)
  library(openxlsx)
  library(flextable)
  library(officer)
  library(broom)
  library(sandwich)
  library(lmtest)
  library(scales)
  library(stringr)
  library(forcats)
})

# -----------------------------------------------------------------------------
# 1. USER SETTINGS
# -----------------------------------------------------------------------------

project_root <- "C:/Users/Stephen/Documents/Researches/Karangura Papers"

paper_dir <- file.path(
  project_root,
  "03_Paper_Water_Governance_Disaster_Resilience"
)

master_csv <- file.path(
  paper_dir,
  "Karangura_BMC_master_cleaned_household_dataset.csv"
)

if (!file.exists(master_csv)) {
  possible_master_files <- list.files(
    paper_dir,
    pattern = "Karangura_BMC_master_cleaned_household_dataset.*\\.csv$|Karangura_BMC_master_cleaned_household_dataset$",
    full.names = TRUE,
    ignore.case = TRUE
  )
  if (length(possible_master_files) > 0) master_csv <- possible_master_files[1]
}

output_dir <- file.path(
  paper_dir,
  paste0("Paper3_Outputs_", format(Sys.Date(), "%Y%m%d"))
)

tables_dir <- file.path(output_dir, "Tables")
figures_dir <- file.path(output_dir, "Figures")
manuscript_dir <- file.path(output_dir, "Manuscript")
objects_dir <- file.path(output_dir, "R_objects")
logs_dir <- file.path(output_dir, "Logs")

purrr::walk(
  c(paper_dir, output_dir, tables_dir, figures_dir, manuscript_dir, objects_dir, logs_dir),
  ~ dir.create(.x, recursive = TRUE, showWarnings = FALSE)
)

if (!file.exists(master_csv)) {
  stop(
    "Master cleaned dataset was not found at:\n", master_csv,
    "\nPlease copy Karangura_BMC_master_cleaned_household_dataset.csv into the Paper 3 folder."
  )
}

# -----------------------------------------------------------------------------
# 2. HELPER FUNCTIONS
# -----------------------------------------------------------------------------

safe_mean <- function(x) if (all(is.na(x))) NA_real_ else mean(x, na.rm = TRUE)
safe_sd <- function(x) if (sum(!is.na(x)) <= 1) NA_real_ else sd(x, na.rm = TRUE)

fmt_n_pct <- function(n, denom, digits = 1) {
  ifelse(is.na(n) | is.na(denom) | denom == 0, NA_character_, paste0(n, " (", round(100 * n / denom, digits), "%)"))
}
fmt_num <- function(x, digits = 2) ifelse(is.na(x), NA_character_, formatC(x, format = "f", digits = digits, big.mark = ","))
fmt_p <- function(p) case_when(is.na(p) ~ NA_character_, p < 0.001 ~ "<0.001", TRUE ~ formatC(p, format = "f", digits = 3))
binary_label <- function(x) case_when(x == 1 ~ "Yes", x == 0 ~ "No", TRUE ~ NA_character_)

safe_chisq_p <- function(x, y) {
  tab <- table(x, y, useNA = "no")
  if (nrow(tab) < 2 | ncol(tab) < 2) return(NA_real_)
  tryCatch(suppressWarnings(chisq.test(tab)$p.value), error = function(e) NA_real_)
}

make_binary_distribution <- function(data, vars, labels) {
  map2_dfr(vars, labels, function(v, lab) {
    if (!v %in% names(data)) {
      return(tibble(variable = v, indicator = lab, denominator = NA_integer_, yes_n = NA_integer_, yes_percent = NA_real_, display = NA_character_))
    }
    x <- data[[v]]
    denom <- sum(!is.na(x))
    yes_n <- sum(x == 1, na.rm = TRUE)
    tibble(variable = v, indicator = lab, denominator = denom, yes_n = yes_n, yes_percent = round(100 * yes_n / denom, 1), display = fmt_n_pct(yes_n, denom))
  })
}

make_group_table <- function(data, group_var, group_label, outcomes, outcome_labels) {
  map2_dfr(outcomes, outcome_labels, function(outcome, outcome_label) {
    if (!all(c(group_var, outcome) %in% names(data))) return(tibble())
    p <- safe_chisq_p(data[[group_var]], data[[outcome]])
    data %>%
      filter(!is.na(.data[[group_var]]), !is.na(.data[[outcome]])) %>%
      group_by(.data[[group_var]]) %>%
      summarise(denominator = n(), yes_n = sum(.data[[outcome]] == 1, na.rm = TRUE), yes_percent = round(100 * yes_n / denominator, 1), .groups = "drop") %>%
      rename(group_category = 1) %>%
      mutate(
        group_category = as.character(group_category),
        group_category = case_when(
          group_category == "1" ~ "Yes",
          group_category == "0" ~ "No",
          TRUE ~ group_category
        ),
        grouping_domain = group_label,
        outcome = outcome_label,
        display = fmt_n_pct(yes_n, denominator),
        p_value = p,
        p_value_display = fmt_p(p)
      ) %>%
      select(grouping_domain, group_category, outcome, denominator, yes_n, yes_percent, display, p_value, p_value_display)
  })
}

modified_poisson <- function(data, outcome, exposures, covariates = NULL, model_label = "Model") {
  all_vars <- c(outcome, exposures, covariates)
  all_vars <- all_vars[all_vars %in% names(data)]
  model_data <- data %>% select(all_of(all_vars)) %>% drop_na()
  if (nrow(model_data) < 30) {
    return(tibble(model = model_label, outcome = outcome, term = NA_character_, prevalence_ratio = NA_real_, conf_low = NA_real_, conf_high = NA_real_, p_value = NA_real_, p_value_display = NA_character_, pr_ci = NA_character_, complete_case_n = nrow(model_data), note = "Model not run: fewer than 30 complete observations"))
  }
  if (length(unique(model_data[[outcome]])) < 2) {
    return(tibble(model = model_label, outcome = outcome, term = NA_character_, prevalence_ratio = NA_real_, conf_low = NA_real_, conf_high = NA_real_, p_value = NA_real_, p_value_display = NA_character_, pr_ci = NA_character_, complete_case_n = nrow(model_data), note = "Model not run: outcome has no variation"))
  }
  rhs <- paste(c(exposures, covariates), collapse = " + ")
  f <- as.formula(paste(outcome, "~", rhs))
  fit <- tryCatch(glm(f, family = poisson(link = "log"), data = model_data), error = function(e) e)
  if (inherits(fit, "error")) {
    return(tibble(model = model_label, outcome = outcome, term = NA_character_, prevalence_ratio = NA_real_, conf_low = NA_real_, conf_high = NA_real_, p_value = NA_real_, p_value_display = NA_character_, pr_ci = NA_character_, complete_case_n = nrow(model_data), note = paste("Model failed:", fit$message)))
  }
  robust <- lmtest::coeftest(fit, vcov. = sandwich::vcovHC(fit, type = "HC0"))
  broom::tidy(robust) %>%
    filter(term != "(Intercept)") %>%
    mutate(model = model_label, outcome = outcome, prevalence_ratio = exp(estimate), conf_low = exp(estimate - 1.96 * std.error), conf_high = exp(estimate + 1.96 * std.error), p_value = p.value, p_value_display = fmt_p(p_value), pr_ci = paste0(fmt_num(prevalence_ratio, 2), " (", fmt_num(conf_low, 2), "-", fmt_num(conf_high, 2), ")"), complete_case_n = nrow(model_data), note = "Modified Poisson regression with robust standard errors") %>%
    select(model, outcome, term, prevalence_ratio, conf_low, conf_high, p_value, p_value_display, pr_ci, complete_case_n, note)
}

# -----------------------------------------------------------------------------
# 3. LOAD AND PREPARE PAPER 3 DATASET
# -----------------------------------------------------------------------------

master <- read_csv(master_csv, show_col_types = FALSE) %>% clean_names()

paper3_data <- master %>%
  filter(paper3_governance_disaster_eligible == 1) %>%
  mutate(
    parish = as.factor(parish),
    respondent_age_group = as.factor(respondent_age_group),
    asset_tertile_label = factor(asset_tertile_label, levels = c("Lowest asset tertile", "Middle asset tertile", "Highest asset tertile")),
    water_governance_enabling = case_when(improved_water_management == 1 & disaster_committee_present == 1 ~ 1, !is.na(improved_water_management) & !is.na(disaster_committee_present) ~ 0, TRUE ~ NA_real_),
    environmental_stewardship_present = case_when(cares_for_environment == 1 | heard_climate_smart_agriculture == 1 ~ 1, cares_for_environment == 0 & heard_climate_smart_agriculture == 0 ~ 0, TRUE ~ NA_real_),
    high_disaster_environment_score = case_when(disaster_environment_score >= 3 ~ 1, disaster_environment_score < 3 ~ 0, TRUE ~ NA_real_),
    disaster_resilience_gap = case_when(low_disaster_environment_resilience == 1 ~ 1, low_disaster_environment_resilience == 0 ~ 0, TRUE ~ NA_real_),
    water_environment_governance_gap = case_when(improved_water_management == 0 | disaster_committee_present == 0 | environmental_stewardship_present == 0 ~ 1, improved_water_management == 1 & disaster_committee_present == 1 & environmental_stewardship_present == 1 ~ 0, TRUE ~ NA_real_),
    high_wash_and_low_resilience = case_when(high_wash_risk == 1 & low_disaster_environment_resilience == 1 ~ 1, !is.na(high_wash_risk) & !is.na(low_disaster_environment_resilience) ~ 0, TRUE ~ NA_real_),
    affected_by_disaster_label = binary_label(affected_by_disaster),
    disaster_committee_label = binary_label(disaster_committee_present),
    improved_water_management_label = binary_label(improved_water_management),
    cares_for_environment_label = binary_label(cares_for_environment),
    low_resilience_label = binary_label(low_disaster_environment_resilience)
  )

analysis_dataset_csv <- file.path(output_dir, "paper3_water_governance_disaster_resilience_analysis_dataset.csv")
write_csv(paper3_data, analysis_dataset_csv)

# -----------------------------------------------------------------------------
# 4. SAMPLE AND DESCRIPTIVE RESULTS
# -----------------------------------------------------------------------------

sample_size <- nrow(paper3_data)

sample_summary <- tibble(
  characteristic = c("Households included in Paper 3 analysis", "Parishes represented", "Villages represented", "Mean respondent age", "Female respondents", "Households in lowest asset tertile", "Vulnerable households, orphanhood or disability", "Households classified as high WASH risk"),
  value = c(as.character(sample_size), as.character(n_distinct(paper3_data$parish, na.rm = TRUE)), as.character(n_distinct(paper3_data$village, na.rm = TRUE)), paste0(fmt_num(safe_mean(paper3_data$respondent_age), 1), " ± ", fmt_num(safe_sd(paper3_data$respondent_age), 1)), fmt_n_pct(sum(paper3_data$respondent_female == 1, na.rm = TRUE), sum(!is.na(paper3_data$respondent_female))), fmt_n_pct(sum(paper3_data$asset_tertile_label == "Lowest asset tertile", na.rm = TRUE), sum(!is.na(paper3_data$asset_tertile_label))), fmt_n_pct(sum(paper3_data$vulnerable_household == 1, na.rm = TRUE), sum(!is.na(paper3_data$vulnerable_household))), fmt_n_pct(sum(paper3_data$high_wash_risk == 1, na.rm = TRUE), sum(!is.na(paper3_data$high_wash_risk))))
)

governance_indicators <- make_binary_distribution(
  paper3_data,
  vars = c("affected_by_disaster", "disaster_committee_present", "improved_water_management", "cares_for_environment", "heard_climate_smart_agriculture", "water_governance_enabling", "environmental_stewardship_present", "high_disaster_environment_score", "low_disaster_environment_resilience", "water_environment_governance_gap", "high_wash_and_low_resilience"),
  labels = c("Household affected by disaster", "Disaster committee present or known", "Improved water management practiced", "Household cares for the environment", "Heard of climate-smart agriculture", "Improved water management plus disaster committee", "Environmental stewardship present", "High disaster/environment score", "Low disaster/environment resilience", "Water/environment governance gap", "High WASH risk plus low resilience")
)

resilience_score_distribution <- paper3_data %>% count(disaster_environment_score, name = "households_n") %>% mutate(percent = round(100 * households_n / sum(households_n), 1), display = fmt_n_pct(households_n, sum(households_n))) %>% arrange(disaster_environment_score)

governance_by_parish <- paper3_data %>%
  group_by(parish) %>%
  summarise(households_n = n(), affected_by_disaster_percent = round(100 * mean(affected_by_disaster == 1, na.rm = TRUE), 1), disaster_committee_percent = round(100 * mean(disaster_committee_present == 1, na.rm = TRUE), 1), improved_water_management_percent = round(100 * mean(improved_water_management == 1, na.rm = TRUE), 1), cares_for_environment_percent = round(100 * mean(cares_for_environment == 1, na.rm = TRUE), 1), high_disaster_environment_score_percent = round(100 * mean(high_disaster_environment_score == 1, na.rm = TRUE), 1), low_resilience_percent = round(100 * mean(low_disaster_environment_resilience == 1, na.rm = TRUE), 1), high_wash_and_low_resilience_percent = round(100 * mean(high_wash_and_low_resilience == 1, na.rm = TRUE), 1), .groups = "drop") %>%
  arrange(desc(low_resilience_percent))

governance_by_asset <- paper3_data %>%
  group_by(asset_tertile_label) %>%
  summarise(households_n = n(), affected_by_disaster_percent = round(100 * mean(affected_by_disaster == 1, na.rm = TRUE), 1), disaster_committee_percent = round(100 * mean(disaster_committee_present == 1, na.rm = TRUE), 1), improved_water_management_percent = round(100 * mean(improved_water_management == 1, na.rm = TRUE), 1), cares_for_environment_percent = round(100 * mean(cares_for_environment == 1, na.rm = TRUE), 1), high_disaster_environment_score_percent = round(100 * mean(high_disaster_environment_score == 1, na.rm = TRUE), 1), low_resilience_percent = round(100 * mean(low_disaster_environment_resilience == 1, na.rm = TRUE), 1), high_wash_and_low_resilience_percent = round(100 * mean(high_wash_and_low_resilience == 1, na.rm = TRUE), 1), .groups = "drop")

outcomes <- c("low_disaster_environment_resilience", "water_environment_governance_gap", "high_wash_and_low_resilience")
outcome_labels <- c("Low disaster/environment resilience", "Water/environment governance gap", "High WASH risk plus low resilience")

group_all <- bind_rows(
  make_group_table(paper3_data, "parish", "Parish", outcomes, outcome_labels),
  make_group_table(paper3_data, "asset_tertile_label", "Household asset tertile", outcomes, outcome_labels),
  make_group_table(paper3_data, "affected_by_disaster_label", "Disaster exposure", outcomes, outcome_labels),
  make_group_table(paper3_data, "vulnerable_household", "Household vulnerability", outcomes, outcome_labels),
  make_group_table(paper3_data, "high_wash_risk", "High WASH risk", outcomes, outcome_labels)
)

# -----------------------------------------------------------------------------
# 5. REGRESSION MODELS
# -----------------------------------------------------------------------------

models_combined <- bind_rows(
  modified_poisson(paper3_data, "low_disaster_environment_resilience", c("high_wash_risk", "affected_by_disaster", "vulnerable_household", "asset_index_simple"), c("respondent_female", "respondent_age", "parish"), "Adjusted model: low disaster/environment resilience"),
  modified_poisson(paper3_data, "water_environment_governance_gap", c("high_wash_risk", "affected_by_disaster", "vulnerable_household", "asset_index_simple"), c("respondent_female", "respondent_age", "parish"), "Adjusted model: water/environment governance gap"),
  modified_poisson(paper3_data, "high_wash_and_low_resilience", c("affected_by_disaster", "vulnerable_household", "asset_index_simple"), c("respondent_female", "respondent_age", "parish"), "Adjusted model: high WASH risk plus low resilience")
)

model_term_labels <- tibble(
  term = c("high_wash_risk", "affected_by_disaster", "vulnerable_household", "asset_index_simple", "respondent_female", "respondent_age"),
  term_label = c("High WASH risk", "Household affected by disaster", "Vulnerable household", "Asset index score", "Female respondent", "Respondent age")
)

models_publication <- models_combined %>%
  left_join(model_term_labels, by = "term") %>%
  mutate(term_label = ifelse(is.na(term_label), term, term_label), `APR (95% CI)` = pr_ci, `p-value` = p_value_display) %>%
  select(Model = model, Outcome = outcome, Variable = term_label, `APR (95% CI)`, `p-value`, `Complete-case N` = complete_case_n, Note = note)

# -----------------------------------------------------------------------------
# 6. PUBLICATION TABLES
# -----------------------------------------------------------------------------

table_1 <- sample_summary
table_2 <- governance_indicators %>% transmute(Indicator = indicator, Denominator = denominator, `Yes, n (%)` = display, `Yes (%)` = yes_percent)
table_3 <- resilience_score_distribution %>% transmute(`Disaster/environment score` = disaster_environment_score, `Households, n (%)` = display, `Households, n` = households_n, Percent = percent)
table_4 <- governance_by_parish
table_5 <- governance_by_asset
table_6 <- group_all %>% transmute(`Grouping domain` = grouping_domain, `Group category` = group_category, Outcome = outcome, Denominator = denominator, `Yes, n (%)` = display, `Yes (%)` = yes_percent, `Chi-square p-value` = p_value_display)
table_7 <- models_publication

table_manuscript_1 <- table_2
table_manuscript_2 <- table_4
table_manuscript_3 <- table_6 %>% filter(`Grouping domain` %in% c("Household asset tertile", "Disaster exposure", "High WASH risk"))
table_manuscript_4 <- table_7


# -----------------------------------------------------------------------------
# 6B. PUBLICATION-READY TABLES FOR JOURNAL SUBMISSION
# -----------------------------------------------------------------------------
# This section keeps the original analysis intact, but adds a publication table
# package. Each table is exported as one Excel sheet and one standalone Word file.

publication_tables_dir <- file.path(output_dir, "Publication_Tables")
publication_word_dir <- file.path(publication_tables_dir, "Standalone_Word_Tables")

dir.create(publication_tables_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(publication_word_dir, recursive = TRUE, showWarnings = FALSE)

format_mean_sd <- function(x) {
  if (all(is.na(x))) return(NA_character_)
  paste0(fmt_num(safe_mean(x), 1), " ± ", fmt_num(safe_sd(x), 1))
}

safe_wilcox_p <- function(x, y) {
  dat <- tibble(x = x, y = y) %>% filter(!is.na(x), !is.na(y))
  if (nrow(dat) < 3 || length(unique(dat$y)) < 2) return(NA_real_)
  tryCatch(suppressWarnings(wilcox.test(x ~ y, data = dat)$p.value), error = function(e) NA_real_)
}

make_pub_continuous_by_outcome <- function(data, vars, labels, outcome, outcome_label) {
  map2_dfr(vars, labels, function(v, lab) {
    if (!all(c(v, outcome) %in% names(data))) return(tibble())
    x <- data[[v]]
    y <- data[[outcome]]
    p <- safe_wilcox_p(x, y)
    tibble(
      Outcome = outcome_label,
      Characteristic = lab,
      Category = "Mean ± SD",
      Overall = format_mean_sd(x),
      `Outcome absent` = format_mean_sd(x[y == 0]),
      `Outcome present` = format_mean_sd(x[y == 1]),
      `p-value` = fmt_p(p),
      `Analysis note` = "Wilcoxon rank-sum test"
    )
  })
}

make_pub_binary_by_outcome <- function(data, vars, labels, outcome, outcome_label) {
  map2_dfr(vars, labels, function(v, lab) {
    if (!all(c(v, outcome) %in% names(data))) return(tibble())
    x <- data[[v]]
    y <- data[[outcome]]
    p <- safe_chisq_p(x, y)

    denom_all <- sum(!is.na(x) & !is.na(y))
    denom_no <- sum(!is.na(x) & y == 0, na.rm = TRUE)
    denom_yes <- sum(!is.na(x) & y == 1, na.rm = TRUE)

    n_all <- sum(x == 1 & !is.na(y), na.rm = TRUE)
    n_no <- sum(x == 1 & y == 0, na.rm = TRUE)
    n_yes <- sum(x == 1 & y == 1, na.rm = TRUE)

    tibble(
      Outcome = outcome_label,
      Characteristic = lab,
      Category = "Yes",
      Overall = fmt_n_pct(n_all, denom_all),
      `Outcome absent` = fmt_n_pct(n_no, denom_no),
      `Outcome present` = fmt_n_pct(n_yes, denom_yes),
      `p-value` = fmt_p(p),
      `Analysis note` = "Chi-square test"
    )
  })
}

make_pub_factor_by_outcome <- function(data, var, label, outcome, outcome_label) {
  if (!all(c(var, outcome) %in% names(data))) return(tibble())

  dat <- data %>%
    select(all_of(c(var, outcome))) %>%
    filter(!is.na(.data[[var]]), !is.na(.data[[outcome]])) %>%
    mutate(group_value = as.character(.data[[var]]))

  if (nrow(dat) == 0) return(tibble())

  p <- safe_chisq_p(dat$group_value, dat[[outcome]])
  levels_found <- sort(unique(dat$group_value))

  map_dfr(levels_found, function(lv) {
    denom_all <- nrow(dat)
    denom_no <- sum(dat[[outcome]] == 0, na.rm = TRUE)
    denom_yes <- sum(dat[[outcome]] == 1, na.rm = TRUE)

    n_all <- sum(dat$group_value == lv, na.rm = TRUE)
    n_no <- sum(dat$group_value == lv & dat[[outcome]] == 0, na.rm = TRUE)
    n_yes <- sum(dat$group_value == lv & dat[[outcome]] == 1, na.rm = TRUE)

    tibble(
      Outcome = outcome_label,
      Characteristic = label,
      Category = lv,
      Overall = fmt_n_pct(n_all, denom_all),
      `Outcome absent` = fmt_n_pct(n_no, denom_no),
      `Outcome present` = fmt_n_pct(n_yes, denom_yes),
      `p-value` = fmt_p(p),
      `Analysis note` = "Chi-square test"
    )
  })
}

make_pub_model_table <- function(model_df, title) {
  model_df %>%
    left_join(model_term_labels, by = "term") %>%
    mutate(
      term_label = ifelse(is.na(term_label), term, term_label),
      Estimate = pr_ci,
      `p-value` = p_value_display,
      Model = model,
      Outcome = case_when(
        outcome == "low_disaster_environment_resilience" ~ "Low disaster/environment resilience",
        outcome == "water_environment_governance_gap" ~ "Water/environment governance gap",
        outcome == "high_wash_and_low_resilience" ~ "High WASH risk plus low resilience",
        TRUE ~ outcome
      )
    ) %>%
    transmute(
      Table = title,
      Model,
      Outcome,
      Variable = term_label,
      `APR (95% CI)` = Estimate,
      `p-value`,
      `Complete-case N` = complete_case_n,
      Note = note
    )
}

pub_table_1_sample_by_low_resilience <- bind_rows(
  make_pub_continuous_by_outcome(
    paper3_data,
    vars = c("respondent_age", "hh_roster_members_n", "hh_under5_members_n", "asset_index_simple"),
    labels = c("Respondent age", "Household members", "Under-five children in household", "Asset index score"),
    outcome = "low_disaster_environment_resilience",
    outcome_label = "Low disaster/environment resilience"
  ),
  make_pub_binary_by_outcome(
    paper3_data,
    vars = c("respondent_female", "orphan_household", "disability_household", "vulnerable_household", "high_wash_risk"),
    labels = c("Female respondent", "Household with orphans", "Household with disability", "Vulnerable household", "High WASH risk"),
    outcome = "low_disaster_environment_resilience",
    outcome_label = "Low disaster/environment resilience"
  ),
  make_pub_factor_by_outcome(paper3_data, "parish", "Parish", "low_disaster_environment_resilience", "Low disaster/environment resilience"),
  make_pub_factor_by_outcome(paper3_data, "asset_tertile_label", "Household asset tertile", "low_disaster_environment_resilience", "Low disaster/environment resilience")
)

pub_table_2_governance_indicators <- table_2 %>%
  transmute(Indicator, Denominator, `Households, n (%)` = `Yes, n (%)`, `Households (%)` = `Yes (%)`)

pub_table_3_score_distribution <- table_3

pub_table_4_parish_profile <- governance_by_parish %>%
  transmute(
    Parish = parish,
    `Households, n` = households_n,
    `Affected by disaster (%)` = affected_by_disaster_percent,
    `Disaster committee present/known (%)` = disaster_committee_percent,
    `Improved water management (%)` = improved_water_management_percent,
    `Cares for environment (%)` = cares_for_environment_percent,
    `High disaster/environment score (%)` = high_disaster_environment_score_percent,
    `Low disaster/environment resilience (%)` = low_resilience_percent,
    `High WASH risk plus low resilience (%)` = high_wash_and_low_resilience_percent
  )

pub_table_5_asset_profile <- governance_by_asset %>%
  transmute(
    `Household asset tertile` = asset_tertile_label,
    `Households, n` = households_n,
    `Affected by disaster (%)` = affected_by_disaster_percent,
    `Disaster committee present/known (%)` = disaster_committee_percent,
    `Improved water management (%)` = improved_water_management_percent,
    `Cares for environment (%)` = cares_for_environment_percent,
    `High disaster/environment score (%)` = high_disaster_environment_score_percent,
    `Low disaster/environment resilience (%)` = low_resilience_percent,
    `High WASH risk plus low resilience (%)` = high_wash_and_low_resilience_percent
  )

pub_table_6_bivariate_associations <- table_6 %>% arrange(`Grouping domain`, Outcome, `Group category`)

pub_table_7_adjusted_models <- make_pub_model_table(models_combined, "Adjusted modified Poisson regression models")

pub_table_8_full_group_matrix <- group_all %>%
  transmute(
    `Grouping domain` = grouping_domain,
    `Group category` = as.character(group_category),
    Outcome = outcome,
    Denominator = denominator,
    `Households, n (%)` = display,
    `Households (%)` = yes_percent,
    `Chi-square p-value` = p_value_display
  ) %>%
  arrange(`Grouping domain`, Outcome, `Group category`)

publication_table_list <- list(
  "Table 1 Sample by low resilience" = pub_table_1_sample_by_low_resilience,
  "Table 2 Governance indicators" = pub_table_2_governance_indicators,
  "Table 3 Score distribution" = pub_table_3_score_distribution,
  "Table 4 Parish profile" = pub_table_4_parish_profile,
  "Table 5 Asset profile" = pub_table_5_asset_profile,
  "Table 6 Bivariate associations" = pub_table_6_bivariate_associations,
  "Table 7 Adjusted models" = pub_table_7_adjusted_models,
  "Table 8 Full group matrix" = pub_table_8_full_group_matrix
)

publication_captions <- c(
  "Table 1 Sample by low resilience" = "Table 1. Baseline household and respondent characteristics by low disaster/environment resilience status",
  "Table 2 Governance indicators" = "Table 2. Water governance, environmental stewardship and disaster-resilience indicators",
  "Table 3 Score distribution" = "Table 3. Distribution of the disaster/environment score",
  "Table 4 Parish profile" = "Table 4. Parish-level governance, environmental stewardship and resilience profile",
  "Table 5 Asset profile" = "Table 5. Governance, environmental stewardship and resilience indicators by household asset tertile",
  "Table 6 Bivariate associations" = "Table 6. Bivariate associations with governance and resilience outcomes",
  "Table 7 Adjusted models" = "Table 7. Adjusted prevalence ratios for governance and resilience outcomes",
  "Table 8 Full group matrix" = "Table 8. Full group comparison matrix for governance and resilience outcomes"
)


# -----------------------------------------------------------------------------
# 7. FIGURES
# -----------------------------------------------------------------------------

theme_pub <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(plot.title = element_text(face = "bold", size = base_size + 3), plot.subtitle = element_text(size = base_size), axis.title = element_text(face = "bold"), axis.text.x = element_text(angle = 30, hjust = 1), panel.grid.minor = element_blank(), legend.position = "bottom", legend.title = element_text(face = "bold"))
}

fig1 <- governance_indicators %>%
  ggplot(aes(x = reorder(indicator, yes_percent), y = yes_percent)) +
  geom_col(width = 0.72, fill = "#1F4E79") +
  geom_text(aes(label = paste0(yes_percent, "%")), hjust = -0.1, size = 3.8) +
  coord_flip() + scale_y_continuous(limits = c(0, 100)) +
  labs(title = "Water governance, environmental stewardship and resilience indicators", subtitle = "Karangura CFCT household baseline", x = NULL, y = "Percentage of households") + theme_pub()
ggsave(file.path(figures_dir, "Figure_1_governance_resilience_indicators.png"), fig1, width = 10.5, height = 7, dpi = 300)

fig2 <- resilience_score_distribution %>%
  ggplot(aes(x = factor(disaster_environment_score), y = percent)) +
  geom_col(fill = "#7030A0", width = 0.7) +
  geom_text(aes(label = paste0(percent, "%")), vjust = -0.35, size = 3.8) +
  labs(title = "Distribution of disaster/environment score", subtitle = "Higher scores indicate stronger enabling conditions", x = "Disaster/environment score", y = "Percentage of households") +
  theme_pub() + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
ggsave(file.path(figures_dir, "Figure_2_disaster_environment_score_distribution.png"), fig2, width = 8, height = 5, dpi = 300)

fig3 <- governance_by_parish %>%
  select(parish, affected_by_disaster_percent, disaster_committee_percent, improved_water_management_percent, cares_for_environment_percent, low_resilience_percent) %>%
  pivot_longer(-parish, names_to = "indicator", values_to = "percent") %>%
  mutate(indicator = case_when(indicator == "affected_by_disaster_percent" ~ "Affected by disaster", indicator == "disaster_committee_percent" ~ "Disaster committee", indicator == "improved_water_management_percent" ~ "Improved water management", indicator == "cares_for_environment_percent" ~ "Cares for environment", indicator == "low_resilience_percent" ~ "Low resilience", TRUE ~ indicator)) %>%
  ggplot(aes(x = parish, y = percent, fill = indicator)) + geom_col(position = position_dodge(width = 0.75), width = 0.68) +
  labs(title = "Governance and resilience indicators by parish", subtitle = "Karangura CFCT household baseline", x = "Parish", y = "Percentage", fill = "Indicator") + theme_pub()
ggsave(file.path(figures_dir, "Figure_3_governance_resilience_by_parish.png"), fig3, width = 10, height = 6, dpi = 300)

fig4 <- governance_by_asset %>%
  select(asset_tertile_label, disaster_committee_percent, improved_water_management_percent, cares_for_environment_percent, low_resilience_percent) %>%
  pivot_longer(-asset_tertile_label, names_to = "indicator", values_to = "percent") %>%
  mutate(indicator = case_when(indicator == "disaster_committee_percent" ~ "Disaster committee", indicator == "improved_water_management_percent" ~ "Improved water management", indicator == "cares_for_environment_percent" ~ "Cares for environment", indicator == "low_resilience_percent" ~ "Low resilience", TRUE ~ indicator)) %>%
  ggplot(aes(x = asset_tertile_label, y = percent, fill = indicator)) + geom_col(position = position_dodge(width = 0.75), width = 0.68) +
  labs(title = "Governance and resilience indicators by asset tertile", subtitle = "Karangura CFCT household baseline", x = "Household asset tertile", y = "Percentage", fill = "Indicator") + theme_pub()
ggsave(file.path(figures_dir, "Figure_4_governance_resilience_by_asset_tertile.png"), fig4, width = 10, height = 6, dpi = 300)

fig5 <- models_combined %>%
  filter(!is.na(prevalence_ratio)) %>%
  left_join(model_term_labels, by = "term") %>%
  mutate(term_label = ifelse(is.na(term_label), term, term_label), outcome_label = case_when(outcome == "low_disaster_environment_resilience" ~ "Low resilience", outcome == "water_environment_governance_gap" ~ "Governance gap", outcome == "high_wash_and_low_resilience" ~ "High WASH risk plus low resilience", TRUE ~ outcome)) %>%
  ggplot(aes(x = prevalence_ratio, y = reorder(term_label, prevalence_ratio))) +
  geom_vline(xintercept = 1, linetype = "dashed") + geom_errorbarh(aes(xmin = conf_low, xmax = conf_high), height = 0.2) + geom_point(size = 2.8) +
  facet_wrap(~ outcome_label, scales = "free_y") + scale_x_log10() +
  labs(title = "Adjusted associations with governance and resilience outcomes", subtitle = "Modified Poisson regression with robust standard errors", x = "Adjusted prevalence ratio, log scale", y = NULL) + theme_pub()
ggsave(file.path(figures_dir, "Figure_5_adjusted_resilience_associations.png"), fig5, width = 11, height = 7, dpi = 300)

# -----------------------------------------------------------------------------
# 8. EXPORT EXCEL
# -----------------------------------------------------------------------------

excel_output <- file.path(tables_dir, "Paper3_Water_Governance_Disaster_Resilience_Tables.xlsx")
wb <- createWorkbook()
sheet_list <- list("Sample summary" = table_1, "Governance indicators" = table_2, "Score distribution" = table_3, "By parish" = table_4, "By asset tertile" = table_5, "Bivariate associations" = table_6, "Regression models" = table_7, "Analysis dataset preview" = head(paper3_data, 100))
header_style <- createStyle(textDecoration = "bold", fgFill = "#1F4E79", fontColour = "white", halign = "center", valign = "center", border = "TopBottomLeftRight", wrapText = TRUE)
body_style <- createStyle(border = "TopBottomLeftRight", valign = "center", wrapText = TRUE)
for (sheet_name in names(sheet_list)) {
  addWorksheet(wb, sheet_name)
  writeData(wb, sheet_name, sheet_list[[sheet_name]])
  n_rows <- max(1, nrow(sheet_list[[sheet_name]]) + 1)
  n_cols <- max(1, ncol(sheet_list[[sheet_name]]))
  addStyle(wb, sheet_name, header_style, rows = 1, cols = 1:n_cols, gridExpand = TRUE)
  if (n_rows >= 2) addStyle(wb, sheet_name, body_style, rows = 2:n_rows, cols = 1:n_cols, gridExpand = TRUE)
  setColWidths(wb, sheet_name, cols = 1:n_cols, widths = "auto")
  freezePane(wb, sheet_name, firstRow = TRUE)
}
saveWorkbook(wb, excel_output, overwrite = TRUE)
rm(wb)
gc()


# Export a second workbook containing only publication-ready tables.
publication_excel_output <- file.path(publication_tables_dir, "Paper3_Publication_Ready_Tables.xlsx")

pub_wb <- createWorkbook()

pub_title_style <- createStyle(textDecoration = "bold", fgFill = "#D9EAF7", fontColour = "#1F4E79", halign = "left", valign = "center", wrapText = TRUE, fontSize = 12)
pub_header_style <- createStyle(textDecoration = "bold", fgFill = "#1F4E79", fontColour = "white", halign = "center", valign = "center", border = "TopBottomLeftRight", wrapText = TRUE)
pub_body_style <- createStyle(border = "TopBottomLeftRight", valign = "center", wrapText = TRUE)
pub_alt_style <- createStyle(fgFill = "#F7F7F7", border = "TopBottomLeftRight", valign = "center", wrapText = TRUE)
pub_note_style <- createStyle(textDecoration = "italic", fontColour = "#666666", wrapText = TRUE)

# Excel sheet names must be 31 characters or fewer.
publication_sheet_names <- c(
  "Table 1 Sample by low resilience" = "T1 Sample low resilience",
  "Table 2 Governance indicators" = "T2 Governance indicators",
  "Table 3 Score distribution" = "T3 Score distribution",
  "Table 4 Parish profile" = "T4 Parish profile",
  "Table 5 Asset profile" = "T5 Asset profile",
  "Table 6 Bivariate associations" = "T6 Bivariate assoc",
  "Table 7 Adjusted models" = "T7 Adjusted models",
  "Table 8 Full group matrix" = "T8 Full group matrix"
)

for (table_name in names(publication_table_list)) {
  sheet_name <- publication_sheet_names[[table_name]]
  addWorksheet(pub_wb, sheet_name)

  pub_dat <- publication_table_list[[table_name]]
  caption <- publication_captions[[table_name]]

  writeData(pub_wb, sheet_name, caption, startRow = 1, startCol = 1)
  mergeCells(pub_wb, sheet_name, cols = 1:max(1, ncol(pub_dat)), rows = 1)
  addStyle(pub_wb, sheet_name, pub_title_style, rows = 1, cols = 1:max(1, ncol(pub_dat)), gridExpand = TRUE)

  writeData(pub_wb, sheet_name, pub_dat, startRow = 3, startCol = 1)

  n_rows <- max(1, nrow(pub_dat) + 3)
  n_cols <- max(1, ncol(pub_dat))

  addStyle(pub_wb, sheet_name, pub_header_style, rows = 3, cols = 1:n_cols, gridExpand = TRUE)

  if (n_rows >= 4) {
    body_rows <- 4:n_rows
    addStyle(pub_wb, sheet_name, pub_body_style, rows = body_rows, cols = 1:n_cols, gridExpand = TRUE)
    alt_rows <- body_rows[seq_along(body_rows) %% 2 == 0]
    if (length(alt_rows) > 0) {
      addStyle(pub_wb, sheet_name, pub_alt_style, rows = alt_rows, cols = 1:n_cols, gridExpand = TRUE)
    }
  }

  note_row <- n_rows + 2
  writeData(pub_wb, sheet_name, "Notes: Percentages are based on non-missing denominators. P-values are from chi-square tests for categorical variables and Wilcoxon rank-sum tests for continuous variables. APR estimates are from modified Poisson regression with robust standard errors.", startRow = note_row, startCol = 1)
  mergeCells(pub_wb, sheet_name, cols = 1:n_cols, rows = note_row)
  addStyle(pub_wb, sheet_name, pub_note_style, rows = note_row, cols = 1:n_cols, gridExpand = TRUE)

  setColWidths(pub_wb, sheet_name, cols = 1:n_cols, widths = "auto")
  freezePane(pub_wb, sheet_name, firstActiveRow = 4)
}

saveWorkbook(pub_wb, publication_excel_output, overwrite = TRUE)
rm(pub_wb)
gc()

# Export each publication-ready table to its own standalone Word document.
make_pub_table_for_word <- function(data, caption) {
  flextable(data) %>%
    set_caption(caption) %>%
    theme_booktabs() %>%
    autofit() %>%
    fontsize(size = 8.5, part = "all") %>%
    bold(part = "header") %>%
    align(align = "center", part = "all") %>%
    align(j = 1, align = "left", part = "body") %>%
    valign(valign = "center", part = "all")
}

safe_file_name <- function(x) {
  x %>%
    str_replace_all("[^A-Za-z0-9]+", "_") %>%
    str_replace_all("_+", "_") %>%
    str_replace_all("^_|_$", "")
}

publication_word_outputs <- character(0)

for (sheet_name in names(publication_table_list)) {
  pub_dat <- publication_table_list[[sheet_name]]
  caption <- publication_captions[[sheet_name]]

  table_file <- file.path(publication_word_dir, paste0("Paper3_", safe_file_name(sheet_name), ".docx"))

  table_doc <- read_docx() %>%
    body_add_par(caption, style = "heading 1") %>%
    body_add_par("Community Water Management, Environmental Stewardship and Disaster Resilience in Karangura, Uganda: Evidence from a Rural Household Baseline Survey", style = "Normal") %>%
    body_add_flextable(make_pub_table_for_word(pub_dat, caption)) %>%
    body_add_par("Notes: Percentages are based on non-missing denominators. P-values are from chi-square tests for categorical variables and Wilcoxon rank-sum tests for continuous variables. APR estimates are from modified Poisson regression with robust standard errors.", style = "Normal")

  print(table_doc, target = table_file)
  publication_word_outputs <- c(publication_word_outputs, table_file)
}

publication_word_index <- file.path(publication_word_dir, "Paper3_Standalone_Word_Tables_Index.txt")
writeLines(publication_word_outputs, publication_word_index)


# -----------------------------------------------------------------------------
# 9. WORD MANUSCRIPT
# -----------------------------------------------------------------------------

get_prev <- function(tbl, var_name) tbl %>% filter(variable == var_name) %>% pull(yes_percent)
prev_affected_disaster <- get_prev(governance_indicators, "affected_by_disaster")
prev_disaster_committee <- get_prev(governance_indicators, "disaster_committee_present")
prev_improved_water_management <- get_prev(governance_indicators, "improved_water_management")
prev_cares_environment <- get_prev(governance_indicators, "cares_for_environment")
prev_heard_csa <- get_prev(governance_indicators, "heard_climate_smart_agriculture")
prev_low_resilience <- get_prev(governance_indicators, "low_disaster_environment_resilience")
prev_high_wash_low_res <- get_prev(governance_indicators, "high_wash_and_low_resilience")

abstract_results <- paste0("Among ", sample_size, " households, ", prev_affected_disaster, "% reported being affected by disaster. Disaster committee presence or awareness was reported by ", prev_disaster_committee, "%, while improved water management was reported by ", prev_improved_water_management, "%. Environmental care was reported by ", prev_cares_environment, "% and awareness of climate-smart agriculture by ", prev_heard_csa, "%. Low disaster/environment resilience was observed in ", prev_low_resilience, "% of households, and ", prev_high_wash_low_res, "% had both high WASH risk and low resilience.")
results_paragraph_1 <- paste0("The analysis included ", sample_size, " households. Disaster exposure was reported by ", prev_affected_disaster, "% of households. Governance-related enabling conditions were uneven: ", prev_disaster_committee, "% reported disaster committee presence or awareness, while ", prev_improved_water_management, "% reported improved water management.")
results_paragraph_2 <- paste0("Environmental stewardship indicators were also mixed. Environmental care was reported by ", prev_cares_environment, "% of households, and climate-smart agriculture awareness by ", prev_heard_csa, "%. Overall, ", prev_low_resilience, "% of households were classified as having low disaster/environment resilience.")
results_paragraph_3 <- paste0("The intersection of WASH risk and resilience was important. ", prev_high_wash_low_res, "% of households had both high WASH risk and low disaster/environment resilience, indicating a group of households that may require integrated WASH, disaster-risk and environmental-health interventions.")
discussion_core <- paste0("This baseline analysis shows that water management, environmental stewardship and disaster resilience are interconnected household-level concerns in Karangura. The findings suggest that WASH risk should not be addressed in isolation from disaster-risk structures and environmental practices. Households with high WASH risk and low resilience represent a priority group for community-based planning, safe water governance and sustainable environmental-health interventions.")

make_pub_table <- function(data, caption) {
  flextable(data) %>% set_caption(caption) %>% theme_booktabs() %>% autofit() %>% fontsize(size = 8.2, part = "all") %>% bold(part = "header") %>% align(align = "center", part = "all") %>% align(j = 1, align = "left", part = "body") %>% valign(valign = "center", part = "all")
}
add_image_if_exists <- function(doc, img_path, title, width = 6.5, height = 4.5) {
  if (file.exists(img_path)) doc %>% body_add_par(title, style = "heading 2") %>% body_add_img(src = img_path, width = width, height = height)
  else doc %>% body_add_par(title, style = "heading 2") %>% body_add_par(paste("Figure not found:", img_path), style = "Normal")
}

ft1 <- make_pub_table(table_manuscript_1, "Table 1. Water governance, environmental stewardship and resilience indicators")
ft2 <- make_pub_table(table_manuscript_2, "Table 2. Governance and resilience indicators by parish")
ft3 <- make_pub_table(table_manuscript_3, "Table 3. Selected bivariate associations with resilience and governance outcomes")
ft4 <- make_pub_table(table_manuscript_4, "Table 4. Adjusted associations with resilience and governance outcomes")

manuscript_output <- file.path(manuscript_dir, "Paper3_Water_Governance_Disaster_Resilience_Manuscript_Draft.docx")
doc <- read_docx()
doc <- doc %>%
  body_add_par("Community Water Management, Environmental Stewardship and Disaster Resilience in Karangura, Uganda: Evidence from a Rural Household Baseline Survey", style = "heading 1") %>%
  body_add_par("Draft manuscript prepared from the Karangura CFCT household baseline survey", style = "Normal") %>%
  body_add_par("Abstract", style = "heading 1") %>%
  body_add_par("Background", style = "heading 2") %>%
  body_add_par("Safe water governance is closely linked with environmental stewardship and disaster resilience, particularly in rural communities exposed to environmental and livelihood shocks. Household baseline data can help identify gaps in community water management and resilience planning.", style = "Normal") %>%
  body_add_par("Methods", style = "heading 2") %>%
  body_add_par("We analysed cross-sectional household baseline survey data from Karangura Sub-County, Uganda. Indicators included disaster exposure, disaster committee presence or awareness, improved water management, environmental care, climate-smart agriculture awareness, disaster/environment resilience and combined high WASH risk plus low resilience. Descriptive statistics, parish and asset comparisons, bivariate tests and modified Poisson regression with robust standard errors were used.", style = "Normal") %>%
  body_add_par("Results", style = "heading 2") %>%
  body_add_par(abstract_results, style = "Normal") %>%
  body_add_par("Conclusion", style = "heading 2") %>%
  body_add_par("The findings show that WASH risk, water-management practices, environmental stewardship and disaster resilience should be considered together in community-level planning. Households with both high WASH risk and low resilience require priority attention.", style = "Normal") %>%
  body_add_par("Introduction", style = "heading 1") %>%
  body_add_par("Water governance, sanitation and environmental management are closely linked in rural public-health settings. Community water management affects household water safety, while disaster exposure and environmental degradation may undermine water access, sanitation and health. Evidence on these household-level relationships is important for integrated planning.", style = "Normal") %>%
  body_add_par("Karangura Sub-County is a rural setting where the Child Focused Community Transformation baseline survey collected household data on water management, disaster exposure, environmental care and related household conditions. This paper examines community water management, environmental stewardship and disaster resilience using household baseline data.", style = "Normal") %>%
  body_add_par("Methods", style = "heading 1") %>%
  body_add_par("Study design and setting", style = "heading 2") %>%
  body_add_par("This was a cross-sectional analysis of household baseline survey data collected in Karangura Sub-County, Kabarole District, Uganda.", style = "Normal") %>%
  body_add_par("Study population", style = "heading 2") %>%
  body_add_par("The analysis included households with available information on improved water management, disaster committee presence or related disaster/environment indicators.", style = "Normal") %>%
  body_add_par("Variables", style = "heading 2") %>%
  body_add_par("Key indicators included disaster exposure, disaster committee presence or awareness, improved water management, care for the environment, climate-smart agriculture awareness and a composite disaster/environment score. Outcomes included low disaster/environment resilience, water/environment governance gap and combined high WASH risk plus low resilience.", style = "Normal") %>%
  body_add_par("Statistical analysis", style = "heading 2") %>%
  body_add_par("Descriptive statistics were used to summarise governance, environmental stewardship and resilience indicators. Comparisons were made by parish, asset tertile, disaster exposure and WASH risk. Modified Poisson regression with robust standard errors was used to estimate adjusted prevalence ratios.", style = "Normal") %>%
  body_add_par("Results", style = "heading 1") %>%
  body_add_par("Governance and resilience indicators", style = "heading 2") %>%
  body_add_par(results_paragraph_1, style = "Normal") %>%
  body_add_par(results_paragraph_2, style = "Normal") %>%
  body_add_par("Intersection of WASH risk and resilience", style = "heading 2") %>%
  body_add_par(results_paragraph_3, style = "Normal") %>%
  body_add_par("Tables", style = "heading 2") %>%
  body_add_flextable(ft1) %>% body_add_par("", style = "Normal") %>%
  body_add_flextable(ft2) %>% body_add_par("", style = "Normal") %>%
  body_add_flextable(ft3) %>% body_add_par("", style = "Normal") %>%
  body_add_flextable(ft4) %>%
  body_add_par("Figures", style = "heading 1")

doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_1_governance_resilience_indicators.png"), "Figure 1. Water governance, environmental stewardship and resilience indicators", 6.5, 4.8)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_2_disaster_environment_score_distribution.png"), "Figure 2. Distribution of disaster/environment score", 6.2, 3.8)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_3_governance_resilience_by_parish.png"), "Figure 3. Governance and resilience indicators by parish", 6.5, 4.2)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_4_governance_resilience_by_asset_tertile.png"), "Figure 4. Governance and resilience indicators by asset tertile", 6.5, 4.2)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_5_adjusted_resilience_associations.png"), "Figure 5. Adjusted associations with resilience and governance outcomes", 6.5, 4.8)

doc <- doc %>%
  body_add_par("Discussion", style = "heading 1") %>%
  body_add_par(discussion_core, style = "Normal") %>%
  body_add_par("The study should be interpreted as a cross-sectional baseline analysis. It identifies household-level patterns and associations but does not establish causality. Some indicators are self-reported and may be affected by recall or reporting bias. Nevertheless, the findings provide useful evidence for integrating WASH, disaster-risk reduction and environmental stewardship in Karangura.", style = "Normal") %>%
  body_add_par("Conclusion", style = "heading 1") %>%
  body_add_par("Water management, environmental stewardship and disaster resilience are important and connected public-health issues in Karangura. Integrated community planning should prioritize households with high WASH risk, weak governance support and low disaster/environment resilience.", style = "Normal") %>%
  body_add_par("Declarations", style = "heading 1") %>%
  body_add_par("Ethics approval and consent to participate: To be completed using the approved Karangura CFCT baseline documentation.", style = "Normal") %>%
  body_add_par("Consent for publication: Not applicable.", style = "Normal") %>%
  body_add_par("Availability of data and materials: To be completed after author and institutional decision on data sharing.", style = "Normal") %>%
  body_add_par("Competing interests: The authors declare no competing interests.", style = "Normal") %>%
  body_add_par("Funding: To be completed.", style = "Normal") %>%
  body_add_par("Authors' contributions: To be completed after authorship confirmation.", style = "Normal") %>%
  body_add_par("Acknowledgements: To be completed.", style = "Normal")

print(doc, target = manuscript_output)

# -----------------------------------------------------------------------------
# 10. SAVE OBJECTS AND LOG
# -----------------------------------------------------------------------------

analysis_objects <- list(
  paper3_data = paper3_data,
  sample_summary = sample_summary,
  governance_indicators = governance_indicators,
  resilience_score_distribution = resilience_score_distribution,
  governance_by_parish = governance_by_parish,
  governance_by_asset = governance_by_asset,
  group_all = group_all,
  models_combined = models_combined,
  models_publication = models_publication,
  publication_table_list = publication_table_list,
  publication_word_outputs = publication_word_outputs,
  tables = list(table_1 = table_1, table_2 = table_2, table_3 = table_3, table_4 = table_4, table_5 = table_5, table_6 = table_6, table_7 = table_7)
)
saveRDS(analysis_objects, file.path(objects_dir, "Paper3_Water_Governance_Disaster_Resilience_Analysis_Objects.rds"))
log_text <- c("PAPER 3 ANALYSIS COMPLETED", paste0("Date/time: ", Sys.time()), paste0("Sample size: ", sample_size), paste0("Output folder: ", output_dir), paste0("Analysis dataset: ", analysis_dataset_csv), paste0("Excel workbook: ", excel_output), paste0("Publication-ready Excel tables: ", publication_excel_output), paste0("Publication-ready standalone Word tables folder: ", publication_word_dir), paste0("Manuscript draft: ", manuscript_output), paste0("Figures directory: ", figures_dir))
writeLines(log_text, file.path(logs_dir, "Paper3_analysis_log.txt"))

cat("\n====================================================================\n")
cat("PAPER 3 ANALYSIS COMPLETED\n")
cat("====================================================================\n\n")
cat("Sample size:", sample_size, "\n")
cat("Output folder:", output_dir, "\n")
cat("Analysis dataset:", analysis_dataset_csv, "\n")
cat("Excel workbook:", excel_output, "\n")
cat("Publication-ready Excel tables:", publication_excel_output, "\n")
cat("Publication-ready standalone Word tables folder:", publication_word_dir, "\n")
cat("Manuscript draft:", manuscript_output, "\n")
cat("Figures:", figures_dir, "\n\n")
cat("Key governance and resilience indicators:\n")
print(table_2)
cat("\nModel results are saved in the Excel workbook and Word draft.\n")
gc()