# =============================================================================
# KARANGURA CFCT BASELINE
# PAPER 5 ANALYSIS SCRIPT
#
# Manuscript:
# Household Energy, Environmental Health and Sanitation Practices in Rural
# Karangura, Uganda:
# Implications for Sustainable Water and Public Health
#
# Input folder:
#   C:/Users/Stephen/Documents/Researches/Karangura Papers/
#     05_Paper_Energy_Environmental_Health
#
# Required input in that folder:
#   Karangura_BMC_master_cleaned_household_dataset.csv
#
# Outputs:
#   Paper5_Outputs_YYYYMMDD/Tables
#   Paper5_Outputs_YYYYMMDD/Figures
#   Paper5_Outputs_YYYYMMDD/Manuscript
#   Paper5_Outputs_YYYYMMDD/R_objects
#   Paper5_Outputs_YYYYMMDD/Logs
#
# Author: Olet Stephen Charles
# =============================================================================


# -----------------------------------------------------------------------------
# 0. CLEAN SESSION AND LOAD PACKAGES
# -----------------------------------------------------------------------------

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,
  "05_Paper_Energy_Environmental_Health"
)

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("Paper5_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 5 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_outcome_by_exposure <- function(data, outcome, outcome_label, exposure, exposure_label) {
  if (!all(c(outcome, exposure) %in% names(data))) return(tibble())

  model_data <- data %>%
    select(all_of(c(outcome, exposure))) %>%
    filter(!is.na(.data[[outcome]]), !is.na(.data[[exposure]])) %>%
    mutate(exposure_group = binary_label(.data[[exposure]]))

  p <- safe_chisq_p(model_data[[exposure]], model_data[[outcome]])

  model_data %>%
    group_by(exposure_group) %>%
    summarise(
      denominator = n(),
      outcome_yes_n = sum(.data[[outcome]] == 1, na.rm = TRUE),
      outcome_percent = round(100 * outcome_yes_n / denominator, 1),
      .groups = "drop"
    ) %>%
    mutate(
      outcome = outcome_label,
      exposure = exposure_label,
      display = fmt_n_pct(outcome_yes_n, denominator),
      p_value = p,
      p_value_display = fmt_p(p)
    ) %>%
    select(
      outcome, exposure, exposure_group,
      denominator, outcome_yes_n, outcome_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 5 DATASET
# -----------------------------------------------------------------------------

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

paper5_data <- master %>%
  filter(paper5_energy_environmental_health_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")
    ),
    poor_sanitation = case_when(
      improved_sanitation == 0 ~ 1,
      improved_sanitation == 1 ~ 0,
      TRUE ~ NA_real_
    ),
    poor_hygiene = case_when(
      handwashing_with_soap_regular == 0 ~ 1,
      handwashing_with_soap_regular == 1 ~ 0,
      TRUE ~ NA_real_
    ),
    household_environmental_health_risk = case_when(
      high_environmental_health_risk == 1 | high_wash_risk == 1 ~ 1,
      high_environmental_health_risk == 0 & high_wash_risk == 0 ~ 0,
      TRUE ~ NA_real_
    ),
    combined_energy_wash_risk = case_when(
      biomass_cooking_fuel == 1 & high_wash_risk == 1 ~ 1,
      !is.na(biomass_cooking_fuel) & !is.na(high_wash_risk) ~ 0,
      TRUE ~ NA_real_
    ),
    combined_energy_sanitation_risk = case_when(
      biomass_cooking_fuel == 1 & poor_sanitation == 1 ~ 1,
      !is.na(biomass_cooking_fuel) & !is.na(poor_sanitation) ~ 0,
      TRUE ~ NA_real_
    ),
    sustainable_household_environment = case_when(
      biomass_cooking_fuel == 0 &
        unsafe_lighting_fuel == 0 &
        improved_sanitation == 1 &
        handwashing_with_soap_regular == 1 &
        high_risk_water_source_any == 0 ~ 1,
      !is.na(biomass_cooking_fuel) &
        !is.na(unsafe_lighting_fuel) &
        !is.na(improved_sanitation) &
        !is.na(handwashing_with_soap_regular) &
        !is.na(high_risk_water_source_any) ~ 0,
      TRUE ~ NA_real_
    )
  )

analysis_dataset_csv <- file.path(output_dir, "paper5_energy_environmental_health_analysis_dataset.csv")
write_csv(paper5_data, analysis_dataset_csv)


# -----------------------------------------------------------------------------
# 4. DESCRIPTIVE RESULTS
# -----------------------------------------------------------------------------

sample_size <- nrow(paper5_data)

sample_summary <- tibble(
  characteristic = c(
    "Households included in Paper 5 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(paper5_data$parish, na.rm = TRUE)),
    as.character(n_distinct(paper5_data$village, na.rm = TRUE)),
    paste0(fmt_num(safe_mean(paper5_data$respondent_age), 1), " ± ", fmt_num(safe_sd(paper5_data$respondent_age), 1)),
    fmt_n_pct(sum(paper5_data$respondent_female == 1, na.rm = TRUE), sum(!is.na(paper5_data$respondent_female))),
    fmt_n_pct(sum(paper5_data$asset_tertile_label == "Lowest asset tertile", na.rm = TRUE), sum(!is.na(paper5_data$asset_tertile_label))),
    fmt_n_pct(sum(paper5_data$vulnerable_household == 1, na.rm = TRUE), sum(!is.na(paper5_data$vulnerable_household))),
    fmt_n_pct(sum(paper5_data$high_wash_risk == 1, na.rm = TRUE), sum(!is.na(paper5_data$high_wash_risk)))
  )
)

environmental_health_indicators <- make_binary_distribution(
  paper5_data,
  vars = c(
    "biomass_cooking_fuel",
    "unsafe_lighting_fuel",
    "high_risk_water_source_any",
    "poor_sanitation",
    "poor_hygiene",
    "high_wash_risk",
    "high_environmental_health_risk",
    "household_environmental_health_risk",
    "combined_energy_wash_risk",
    "combined_energy_sanitation_risk",
    "sustainable_household_environment"
  ),
  labels = c(
    "Biomass cooking fuel",
    "Unsafe lighting fuel",
    "Any high-risk water source",
    "Poor sanitation",
    "Poor hygiene, no regular handwashing with soap",
    "High WASH risk",
    "High environmental-health risk",
    "High environmental-health or WASH risk",
    "Biomass fuel plus high WASH risk",
    "Biomass fuel plus poor sanitation",
    "Sustainable household environment"
  )
)

environmental_risk_score_distribution <- paper5_data %>%
  count(environmental_health_risk_score, name = "households_n") %>%
  mutate(
    percent = round(100 * households_n / sum(households_n), 1),
    display = fmt_n_pct(households_n, sum(households_n))
  ) %>%
  arrange(environmental_health_risk_score)

energy_wash_by_parish <- paper5_data %>%
  group_by(parish) %>%
  summarise(
    households_n = n(),
    biomass_cooking_fuel_percent = round(100 * mean(biomass_cooking_fuel == 1, na.rm = TRUE), 1),
    unsafe_lighting_fuel_percent = round(100 * mean(unsafe_lighting_fuel == 1, na.rm = TRUE), 1),
    high_risk_water_source_percent = round(100 * mean(high_risk_water_source_any == 1, na.rm = TRUE), 1),
    poor_sanitation_percent = round(100 * mean(poor_sanitation == 1, na.rm = TRUE), 1),
    poor_hygiene_percent = round(100 * mean(poor_hygiene == 1, na.rm = TRUE), 1),
    high_environmental_health_risk_percent = round(100 * mean(high_environmental_health_risk == 1, na.rm = TRUE), 1),
    combined_energy_wash_risk_percent = round(100 * mean(combined_energy_wash_risk == 1, na.rm = TRUE), 1),
    sustainable_household_environment_percent = round(100 * mean(sustainable_household_environment == 1, na.rm = TRUE), 1),
    .groups = "drop"
  ) %>%
  arrange(desc(high_environmental_health_risk_percent))

energy_wash_by_asset <- paper5_data %>%
  group_by(asset_tertile_label) %>%
  summarise(
    households_n = n(),
    biomass_cooking_fuel_percent = round(100 * mean(biomass_cooking_fuel == 1, na.rm = TRUE), 1),
    unsafe_lighting_fuel_percent = round(100 * mean(unsafe_lighting_fuel == 1, na.rm = TRUE), 1),
    high_risk_water_source_percent = round(100 * mean(high_risk_water_source_any == 1, na.rm = TRUE), 1),
    poor_sanitation_percent = round(100 * mean(poor_sanitation == 1, na.rm = TRUE), 1),
    poor_hygiene_percent = round(100 * mean(poor_hygiene == 1, na.rm = TRUE), 1),
    high_environmental_health_risk_percent = round(100 * mean(high_environmental_health_risk == 1, na.rm = TRUE), 1),
    combined_energy_wash_risk_percent = round(100 * mean(combined_energy_wash_risk == 1, na.rm = TRUE), 1),
    sustainable_household_environment_percent = round(100 * mean(sustainable_household_environment == 1, na.rm = TRUE), 1),
    .groups = "drop"
  )


# -----------------------------------------------------------------------------
# 5. BIVARIATE ASSOCIATIONS
# -----------------------------------------------------------------------------

outcomes <- c(
  "high_environmental_health_risk",
  "household_environmental_health_risk",
  "combined_energy_wash_risk",
  "sustainable_household_environment"
)

outcome_labels <- c(
  "High environmental-health risk",
  "High environmental-health or WASH risk",
  "Biomass fuel plus high WASH risk",
  "Sustainable household environment"
)

exposures <- c(
  "biomass_cooking_fuel",
  "unsafe_lighting_fuel",
  "poor_sanitation",
  "poor_hygiene",
  "high_risk_water_source_any",
  "vulnerable_household"
)

exposure_labels <- c(
  "Biomass cooking fuel",
  "Unsafe lighting fuel",
  "Poor sanitation",
  "Poor hygiene",
  "Any high-risk water source",
  "Vulnerable household"
)

bivariate_combined <- map2_dfr(outcomes, outcome_labels, function(outcome, outcome_label) {
  map2_dfr(exposures, exposure_labels, function(exposure, exposure_label) {
    make_outcome_by_exposure(
      paper5_data,
      outcome = outcome,
      outcome_label = outcome_label,
      exposure = exposure,
      exposure_label = exposure_label
    )
  })
})


# -----------------------------------------------------------------------------
# 6. MODIFIED POISSON REGRESSION MODELS
# -----------------------------------------------------------------------------

model_high_environmental_health <- modified_poisson(
  data = paper5_data,
  outcome = "high_environmental_health_risk",
  exposures = c(
    "biomass_cooking_fuel",
    "unsafe_lighting_fuel",
    "poor_sanitation",
    "poor_hygiene",
    "high_risk_water_source_any"
  ),
  covariates = c(
    "respondent_female",
    "respondent_age",
    "vulnerable_household",
    "asset_index_simple",
    "parish"
  ),
  model_label = "Adjusted model: high environmental-health risk"
)

model_combined_energy_wash <- modified_poisson(
  data = paper5_data,
  outcome = "combined_energy_wash_risk",
  exposures = c(
    "unsafe_lighting_fuel",
    "poor_sanitation",
    "poor_hygiene",
    "high_risk_water_source_any"
  ),
  covariates = c(
    "respondent_female",
    "respondent_age",
    "vulnerable_household",
    "asset_index_simple",
    "parish"
  ),
  model_label = "Adjusted model: biomass fuel plus high WASH risk"
)

model_sustainable_environment <- modified_poisson(
  data = paper5_data,
  outcome = "sustainable_household_environment",
  exposures = c(
    "vulnerable_household",
    "asset_index_simple",
    "respondent_female"
  ),
  covariates = c(
    "respondent_age",
    "parish"
  ),
  model_label = "Adjusted model: sustainable household environment"
)

models_combined <- bind_rows(
  model_high_environmental_health,
  model_combined_energy_wash,
  model_sustainable_environment
)

model_term_labels <- tibble(
  term = c(
    "biomass_cooking_fuel",
    "unsafe_lighting_fuel",
    "poor_sanitation",
    "poor_hygiene",
    "high_risk_water_source_any",
    "vulnerable_household",
    "asset_index_simple",
    "respondent_female",
    "respondent_age"
  ),
  term_label = c(
    "Biomass cooking fuel",
    "Unsafe lighting fuel",
    "Poor sanitation",
    "Poor hygiene",
    "Any high-risk water source",
    "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
  )


# -----------------------------------------------------------------------------
# 7. PUBLICATION TABLES
# -----------------------------------------------------------------------------

table_1 <- sample_summary

table_2 <- environmental_health_indicators %>%
  transmute(
    Indicator = indicator,
    Denominator = denominator,
    `Yes, n (%)` = display,
    `Yes (%)` = yes_percent
  )

table_3 <- environmental_risk_score_distribution %>%
  transmute(
    `Environmental-health risk score` = environmental_health_risk_score,
    `Households, n (%)` = display,
    `Households, n` = households_n,
    Percent = percent
  )

table_4 <- energy_wash_by_parish

table_5 <- energy_wash_by_asset

table_6 <- bivariate_combined %>%
  transmute(
    Outcome = outcome,
    Exposure = exposure,
    `Exposure category` = exposure_group,
    Denominator = denominator,
    `Outcome, n (%)` = display,
    `Outcome (%)` = outcome_percent,
    `Chi-square p-value` = p_value_display
  )

table_7 <- models_publication

table_manuscript_1 <- table_2
table_manuscript_2 <- table_3
table_manuscript_3 <- table_4
table_manuscript_4 <- table_7


# -----------------------------------------------------------------------------
# 8. 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 <- environmental_health_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.7) +
  coord_flip() +
  scale_y_continuous(limits = c(0, 100)) +
  labs(
    title = "Household energy and environmental-health indicators",
    subtitle = "Karangura CFCT household baseline",
    x = NULL,
    y = "Percentage of households"
  ) +
  theme_pub()

ggsave(file.path(figures_dir, "Figure_1_energy_environmental_health_indicators.png"), fig1, width = 10.5, height = 7, dpi = 300)

fig2 <- environmental_risk_score_distribution %>%
  ggplot(aes(x = factor(environmental_health_risk_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 environmental-health risk score",
    subtitle = "Higher scores represent more household environmental-health risks",
    x = "Environmental-health risk 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_environmental_health_risk_score_distribution.png"), fig2, width = 8, height = 5, dpi = 300)

fig3 <- energy_wash_by_parish %>%
  select(
    parish,
    biomass_cooking_fuel_percent,
    unsafe_lighting_fuel_percent,
    high_risk_water_source_percent,
    high_environmental_health_risk_percent,
    sustainable_household_environment_percent
  ) %>%
  pivot_longer(-parish, names_to = "indicator", values_to = "percent") %>%
  mutate(
    indicator = case_when(
      indicator == "biomass_cooking_fuel_percent" ~ "Biomass cooking fuel",
      indicator == "unsafe_lighting_fuel_percent" ~ "Unsafe lighting fuel",
      indicator == "high_risk_water_source_percent" ~ "High-risk water source",
      indicator == "high_environmental_health_risk_percent" ~ "High environmental-health risk",
      indicator == "sustainable_household_environment_percent" ~ "Sustainable household environment",
      TRUE ~ indicator
    )
  ) %>%
  ggplot(aes(x = parish, y = percent, fill = indicator)) +
  geom_col(position = position_dodge(width = 0.75), width = 0.68) +
  labs(
    title = "Energy and environmental-health indicators by parish",
    subtitle = "Karangura CFCT household baseline",
    x = "Parish",
    y = "Percentage",
    fill = "Indicator"
  ) +
  theme_pub()

ggsave(file.path(figures_dir, "Figure_3_energy_environmental_health_by_parish.png"), fig3, width = 10.5, height = 6.2, dpi = 300)

fig4 <- energy_wash_by_asset %>%
  select(
    asset_tertile_label,
    biomass_cooking_fuel_percent,
    unsafe_lighting_fuel_percent,
    high_environmental_health_risk_percent,
    combined_energy_wash_risk_percent,
    sustainable_household_environment_percent
  ) %>%
  pivot_longer(-asset_tertile_label, names_to = "indicator", values_to = "percent") %>%
  mutate(
    indicator = case_when(
      indicator == "biomass_cooking_fuel_percent" ~ "Biomass cooking fuel",
      indicator == "unsafe_lighting_fuel_percent" ~ "Unsafe lighting fuel",
      indicator == "high_environmental_health_risk_percent" ~ "High environmental-health risk",
      indicator == "combined_energy_wash_risk_percent" ~ "Biomass plus high WASH risk",
      indicator == "sustainable_household_environment_percent" ~ "Sustainable household environment",
      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 = "Energy and environmental-health 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_energy_environmental_health_by_asset_tertile.png"), fig4, width = 10.5, height = 6.2, dpi = 300)

fig5 <- paper5_data %>%
  mutate(
    energy_wash_group = case_when(
      biomass_cooking_fuel == 0 & high_wash_risk == 0 ~ "Non-biomass and low WASH risk",
      biomass_cooking_fuel == 1 & high_wash_risk == 0 ~ "Biomass only",
      biomass_cooking_fuel == 0 & high_wash_risk == 1 ~ "High WASH risk only",
      biomass_cooking_fuel == 1 & high_wash_risk == 1 ~ "Biomass and high WASH risk",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(energy_wash_group)) %>%
  group_by(energy_wash_group) %>%
  summarise(
    households_n = n(),
    household_diarrhea_percent = round(100 * mean(household_diarrhea_recent == 1, na.rm = TRUE), 1),
    high_environmental_health_risk_percent = round(100 * mean(high_environmental_health_risk == 1, na.rm = TRUE), 1),
    sustainable_household_environment_percent = round(100 * mean(sustainable_household_environment == 1, na.rm = TRUE), 1),
    .groups = "drop"
  ) %>%
  pivot_longer(
    cols = c(household_diarrhea_percent, high_environmental_health_risk_percent, sustainable_household_environment_percent),
    names_to = "outcome",
    values_to = "percent"
  ) %>%
  mutate(
    outcome = case_when(
      outcome == "household_diarrhea_percent" ~ "Recent household diarrhoea",
      outcome == "high_environmental_health_risk_percent" ~ "High environmental-health risk",
      outcome == "sustainable_household_environment_percent" ~ "Sustainable household environment",
      TRUE ~ outcome
    )
  ) %>%
  ggplot(aes(x = energy_wash_group, y = percent, fill = outcome)) +
  geom_col(position = position_dodge(width = 0.75), width = 0.68) +
  labs(
    title = "Health and environmental outcomes by combined energy-WASH risk group",
    subtitle = "Karangura CFCT household baseline",
    x = NULL,
    y = "Percentage",
    fill = "Outcome"
  ) +
  theme_pub()

ggsave(file.path(figures_dir, "Figure_5_outcomes_by_combined_energy_wash_group.png"), fig5, width = 11, height = 6.3, dpi = 300)

fig6 <- bivariate_combined %>%
  filter(Outcome %in% c("High environmental-health risk", "Biomass fuel plus high WASH risk", "Sustainable household environment")) %>%
  filter(Exposure %in% c("Biomass cooking fuel", "Unsafe lighting fuel", "Poor sanitation", "Poor hygiene", "Any high-risk water source")) %>%
  ggplot(aes(x = Exposure, y = `Outcome (%)`, fill = `Exposure category`)) +
  geom_col(position = position_dodge(width = 0.75), width = 0.68) +
  facet_wrap(~ Outcome) +
  labs(
    title = "Environmental-health outcomes by exposure category",
    subtitle = "Bivariate comparison of energy, WASH and sanitation exposures",
    x = NULL,
    y = "Outcome percentage",
    fill = "Exposure category"
  ) +
  theme_pub()

ggsave(file.path(figures_dir, "Figure_6_environmental_outcomes_by_exposure.png"), fig6, width = 12, height = 7, dpi = 300)

fig7 <- 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 == "high_environmental_health_risk" ~ "High environmental-health risk",
      outcome == "combined_energy_wash_risk" ~ "Biomass plus high WASH risk",
      outcome == "sustainable_household_environment" ~ "Sustainable household environment",
      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 environmental-health 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_7_adjusted_environmental_health_associations.png"), fig7, width = 11, height = 7, dpi = 300)


# -----------------------------------------------------------------------------
# 9. EXPORT EXCEL
# -----------------------------------------------------------------------------

excel_output <- file.path(tables_dir, "Paper5_Energy_Environmental_Health_Tables.xlsx")

wb <- createWorkbook()

sheet_list <- list(
  "Sample summary" = table_1,
  "Environmental indicators" = table_2,
  "Risk 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(paper5_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()


# -----------------------------------------------------------------------------
# 10. MANUSCRIPT NARRATIVE
# -----------------------------------------------------------------------------

get_prev <- function(tbl, var_name) {
  tbl %>% filter(variable == var_name) %>% pull(yes_percent)
}

prev_biomass <- get_prev(environmental_health_indicators, "biomass_cooking_fuel")
prev_unsafe_lighting <- get_prev(environmental_health_indicators, "unsafe_lighting_fuel")
prev_high_risk_water <- get_prev(environmental_health_indicators, "high_risk_water_source_any")
prev_poor_sanitation <- get_prev(environmental_health_indicators, "poor_sanitation")
prev_poor_hygiene <- get_prev(environmental_health_indicators, "poor_hygiene")
prev_high_wash <- get_prev(environmental_health_indicators, "high_wash_risk")
prev_high_env <- get_prev(environmental_health_indicators, "high_environmental_health_risk")
prev_combined_energy_wash <- get_prev(environmental_health_indicators, "combined_energy_wash_risk")
prev_sustainable <- get_prev(environmental_health_indicators, "sustainable_household_environment")

abstract_results <- paste0(
  "Among ", sample_size, " households, biomass cooking fuel was reported by ",
  prev_biomass, "% and unsafe lighting fuel by ", prev_unsafe_lighting,
  "%. High-risk water source use was reported by ", prev_high_risk_water,
  "%, poor sanitation by ", prev_poor_sanitation,
  "% and poor hygiene by ", prev_poor_hygiene,
  "%. Overall, ", prev_high_env,
  "% were classified as high environmental-health risk, while ",
  prev_combined_energy_wash, "% had both biomass cooking fuel and high WASH risk. Only ",
  prev_sustainable, "% met the sustainable household environment indicator."
)

results_paragraph_1 <- paste0(
  "The analysis included ", sample_size, " households. Biomass cooking fuel was reported by ",
  prev_biomass, "% of households and unsafe lighting fuel by ",
  prev_unsafe_lighting, "%. These energy indicators coexisted with WASH-related risks, including high-risk water source use in ",
  prev_high_risk_water, "%, poor sanitation in ", prev_poor_sanitation,
  "% and poor hygiene in ", prev_poor_hygiene, "% of households."
)

results_paragraph_2 <- paste0(
  "Overall, ", prev_high_env,
  "% of households were classified as having high environmental-health risk. Combined biomass cooking fuel and high WASH risk was observed in ",
  prev_combined_energy_wash, "% of households. Only ", prev_sustainable,
  "% met the sustainable household environment indicator, which combined cleaner energy, safer water, sanitation and hygiene conditions."
)

discussion_core <- paste0(
  "This baseline analysis shows that household energy, sanitation, hygiene and water-safety conditions are closely linked environmental-health issues in Karangura. ",
  "The findings suggest that sustainable household health cannot be addressed through WASH interventions alone. ",
  "Integrated programming should consider cooking fuel, lighting fuel, water safety, sanitation and handwashing as part of a broader rural environmental-health agenda."
)


# -----------------------------------------------------------------------------
# 11. WORD MANUSCRIPT
# -----------------------------------------------------------------------------

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 <- doc %>%
      body_add_par(title, style = "heading 2") %>%
      body_add_img(src = img_path, width = width, height = height)
  } else {
    doc <- doc %>%
      body_add_par(title, style = "heading 2") %>%
      body_add_par(paste("Figure not found:", img_path), style = "Normal")
  }

  doc
}

ft1 <- make_pub_table(table_manuscript_1, "Table 1. Household energy and environmental-health indicators")
ft2 <- make_pub_table(table_manuscript_2, "Table 2. Environmental-health risk score distribution")
ft3 <- make_pub_table(table_manuscript_3, "Table 3. Energy and environmental-health indicators by parish")
ft4 <- make_pub_table(table_manuscript_4, "Table 4. Adjusted associations with environmental-health outcomes")

manuscript_output <- file.path(
  manuscript_dir,
  "Paper5_Energy_Environmental_Health_Manuscript_Draft.docx"
)

doc <- read_docx()

doc <- doc %>%
  body_add_par("Household Energy, Environmental Health and Sanitation Practices in Rural Karangura, Uganda: Implications for Sustainable Water and Public Health", 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(
    "Household environmental health is shaped by energy use, water safety, sanitation and hygiene. In rural settings, these risks may cluster within the same households, increasing the need for integrated public-health and sustainable-development responses.",
    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. Key indicators included biomass cooking fuel, unsafe lighting fuel, high-risk water source use, poor sanitation, poor hygiene, high WASH risk, high environmental-health risk and a sustainable household environment indicator. Descriptive statistics, parish and asset comparisons, bivariate analysis 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(
    "Environmental-health risks in Karangura extend beyond water and sanitation alone. Household energy, water safety, sanitation and hygiene should be addressed together in sustainable public-health programming.",
    style = "Normal"
  ) %>%

  body_add_par("Introduction", style = "heading 1") %>%
  body_add_par(
    "Household environmental health is influenced by the combined conditions in which families cook, light their homes, obtain water, practise sanitation and maintain hygiene. In rural settings, dependence on biomass cooking fuel, unsafe lighting sources, unsafe water and poor sanitation may cluster within the same households.",
    style = "Normal"
  ) %>%
  body_add_par(
    "Karangura Sub-County is a rural setting where the Child Focused Community Transformation baseline survey collected household information on energy, water, sanitation, hygiene and environmental conditions. This paper examines household energy, environmental-health and sanitation practices and their implications for sustainable water and public health.",
    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 household energy use, sanitation, hygiene and water-safety indicators.",
    style = "Normal"
  ) %>%
  body_add_par("Variables", style = "heading 2") %>%
  body_add_par(
    "Key indicators included biomass cooking fuel, unsafe lighting fuel, high-risk water source use, poor sanitation, poor hygiene, high WASH risk, high environmental-health risk, combined energy-WASH risk and sustainable household environment.",
    style = "Normal"
  ) %>%
  body_add_par("Statistical analysis", style = "heading 2") %>%
  body_add_par(
    "Descriptive statistics were used to summarise household environmental-health indicators. Indicators were compared by parish and household asset tertile. Bivariate associations were assessed using cross-tabulations and chi-square tests. 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("Household energy and environmental-health risks", style = "heading 2") %>%
  body_add_par(results_paragraph_1, style = "Normal") %>%
  body_add_par("Combined environmental-health burden", style = "heading 2") %>%
  body_add_par(results_paragraph_2, 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_energy_environmental_health_indicators.png"), "Figure 1. Household energy and environmental-health indicators", 6.5, 4.8)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_2_environmental_health_risk_score_distribution.png"), "Figure 2. Environmental-health risk score distribution", 6.2, 3.8)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_3_energy_environmental_health_by_parish.png"), "Figure 3. Energy and environmental-health indicators by parish", 6.5, 4.2)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_4_energy_environmental_health_by_asset_tertile.png"), "Figure 4. Energy and environmental-health indicators by asset tertile", 6.5, 4.2)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_5_outcomes_by_combined_energy_wash_group.png"), "Figure 5. Outcomes by combined energy-WASH risk group", 6.5, 4.2)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_6_environmental_outcomes_by_exposure.png"), "Figure 6. Environmental-health outcomes by exposure category", 6.5, 4.4)
doc <- add_image_if_exists(doc, file.path(figures_dir, "Figure_7_adjusted_environmental_health_associations.png"), "Figure 7. Adjusted associations with environmental-health 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 patterns and associations but does not establish causality. Some indicators were self-reported and may be affected by recall or reporting bias. Nevertheless, the results provide practical evidence for integrating household energy, WASH and environmental-health programming.",
    style = "Normal"
  ) %>%
  body_add_par("Conclusion", style = "heading 1") %>%
  body_add_par(
    "Household energy, water safety, sanitation and hygiene risks cluster in Karangura households. Sustainable public-health programming should integrate cleaner household energy, safe water access, sanitation improvement and hygiene promotion.",
    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)


# -----------------------------------------------------------------------------
# 12. SAVE OBJECTS AND LOG
# -----------------------------------------------------------------------------

analysis_objects <- list(
  paper5_data = paper5_data,
  sample_summary = sample_summary,
  environmental_health_indicators = environmental_health_indicators,
  environmental_risk_score_distribution = environmental_risk_score_distribution,
  energy_wash_by_parish = energy_wash_by_parish,
  energy_wash_by_asset = energy_wash_by_asset,
  bivariate_combined = bivariate_combined,
  models_combined = models_combined,
  models_publication = models_publication,
  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, "Paper5_Energy_Environmental_Health_Analysis_Objects.rds")
)

log_text <- c(
  "PAPER 5 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("Manuscript draft: ", manuscript_output),
  paste0("Figures directory: ", figures_dir)
)

writeLines(log_text, file.path(logs_dir, "Paper5_analysis_log.txt"))


# -----------------------------------------------------------------------------
# 13. CONSOLE SUMMARY
# -----------------------------------------------------------------------------

cat("\n====================================================================\n")
cat("PAPER 5 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("Manuscript draft:", manuscript_output, "\n")
cat("Figures:", figures_dir, "\n\n")

cat("Environmental-health indicators:\n")
print(table_2)

cat("\nModel results are saved in the Excel workbook and Word draft.\n")

gc()