# ==============================================================================
# HawaSpatial Pro: Holistic & Areal Weighted Analysis Platform
# Final Publication Version: Fixed Multilevel Engine & 6-Outcome Support
# ==============================================================================
library(shiny)
library(tidyverse)
library(sf)
library(spdep)
library(haven)
library(performance)
library(broom.mixed)
library(bs4Dash)
library(RColorBrewer)
library(gstat)
library(spatialreg)
library(viridis)
library(biscale)
library(scales)
library(ggrepel)
library(glmmTMB)
library(rmarkdown)
library(cowplot)
library(classInt)
library(corrplot)
library(survey)
library(sae)
library(spatialreg)
library(lmtest)
library(ineq)
# Disable S2 to prevent topology errors
sf_use_s2(FALSE)
options(survey.lonely.psu = "adjust")
options(shiny.maxRequestSize = 900 * 1024^2)

custom_theme <- "
  .main-header .logo { background-color: #006400 !important; color: #FFD700 !important; font-weight: bold; }
  .main-header .navbar { background-color: #006400 !important; }
  .sidebar-dark-primary { background-color: #004d00 !important; }
  .btn-success { background-color: #006400 !important; border-color: #FFD700 !important; color: #FFD700 !important; }
  .btn-info { background-color: #004d00 !important; border-color: #FFD700 !important; }
  .nav-pills .nav-link.active, .nav-pills .show > .nav-link { background-color: #006400 !important; color: #FFD700 !important; }
"

# --- Helper Function: Robust ICC for glmmTMB (Used by Shiny UI) ---
calc_freq_icc <- function(model) {
  tryCatch({
    # 1. Extract Cluster Variance (Tau^2)
    # glmmTMB stores VarCorr in a nested list structure $cond$cluster_idx
    vc <- VarCorr(model)
    var_u <- as.numeric(vc$cond$cluster_idx[1,1]) # Extract scalar variance

    # 2. Determine Residual Variance (Sigma^2) based on Family
    fam <- family(model)$family

    if (grepl("binomial", fam)) {
      var_e <- 3.289868
    } else if (grepl("poisson", fam) || grepl("nbinom", fam)) {
      var_e <- 1.0
    } else if (grepl("beta", fam)) {
      var_e <- 1.0
    } else {
      var_e <- sigma(model)^2
    }

    # 3. Calculate ICC
    icc_val <- var_u / (var_u + var_e)

    if (is.na(icc_val) || is.null(icc_val)) return(0)
    return(icc_val)

  }, error = function(e) {
    # If anything fails, return 0 (but print error to console for debugging)
    print(paste("ICC Error:", e$message))
    return(0)
  })
}

# --- Helper Function: Calculate PCV (Proportional Change in Variance) ---
calc_pcv <- function(m0, m_full) {
  icc_m0 <- calc_freq_icc(m0)
  icc_m_full <- calc_freq_icc(m_full)
  if (icc_m0 > 0) {
    pcv_val <- (icc_m0 - icc_m_full) / icc_m0
    return(max(0, pcv_val))
  } else { return(0) }
}



calc_descriptive_table2 <- function(data1, data2, label1, label2, preds, wgt_var) {
  res_list <- list()
  for(p in preds) {
    # Group 1
    df1 <- data1 %>% group_by(Category = as.character(haven::as_factor(!!sym(p)))) %>%
      summarise(n = n(), w_n = sum(!!sym(wgt_var), na.rm=TRUE)) %>%
      mutate(Pct = (w_n/sum(w_n))*100, Group = label1, Variable = p)
    # Group 2
    df2 <- data2 %>% group_by(Category = as.character(haven::as_factor(!!sym(p)))) %>%
      summarise(n = n(), w_n = sum(!!sym(wgt_var), na.rm=TRUE)) %>%
      mutate(Pct = (w_n/sum(w_n))*100, Group = label2, Variable = p)

    combined <- full_join(df1, df2, by=c("Variable", "Category"), suffix=c(paste0("_", label1), paste0("_", label2)))
    res_list[[p]] <- combined
  }
  do.call(rbind, res_list) %>% dplyr::select(Variable, Category, contains("Pct"))
}


# --- ADDED FOR FINAL REPUBLISH VERSION: APA helpers and bivariate screening ---
format_p_apa <- function(p) {
  if (is.null(p) || length(p) == 0 || is.na(p)) return("-")
  if (p < 0.001) return("< .001")
  sub("^0", "", sprintf("%.3f", p))
}

sig_label <- function(p) {
  if (is.null(p) || length(p) == 0 || is.na(p)) return("-")
  if (p < 0.001) return("***")
  if (p < 0.01) return("**")
  if (p < 0.05) return("*")
  "ns"
}

interpret_p <- function(p) {
  if (is.null(p) || length(p) == 0 || is.na(p)) return("Not estimated")
  if (p < 0.05) return("Statistically significant")
  "Not statistically significant"
}

safe_numeric <- function(x) {
  suppressWarnings(as.numeric(as.character(x)))
}

calc_bivariate_screening_table <- function(df, dv_col, cont_vars = NULL, cat_vars = NULL) {
  res <- list()
  if (is.null(df) || !(dv_col %in% names(df))) return(data.frame(Message = "Outcome variable not available."))
  all_cat <- unique(cat_vars[!is.na(cat_vars) & cat_vars %in% names(df)])
  all_cont <- unique(cont_vars[!is.na(cont_vars) & cont_vars %in% names(df)])
  y <- df[[dv_col]]
  y_fac <- as.factor(y)

  for (v in all_cat) {
    tmp <- data.frame(y = y_fac, x = haven::as_factor(df[[v]])) %>% tidyr::drop_na()
    if (nrow(tmp) < 3 || dplyr::n_distinct(tmp$y) < 2 || dplyr::n_distinct(tmp$x) < 2) {
      res[[length(res) + 1]] <- data.frame(
        Variable = v, Variable_Type = "Categorical / factor", Test = "Pearson chi-square",
        Statistic = "Not estimated", df = "-", P_value = "-", Significance = "-",
        Note = "Insufficient variation or sample size"
      )
      next
    }
    tab <- table(tmp$y, tmp$x)
    tst <- tryCatch(suppressWarnings(stats::chisq.test(tab)), error = function(e) NULL)
    if (is.null(tst)) {
      res[[length(res) + 1]] <- data.frame(
        Variable = v, Variable_Type = "Categorical / factor", Test = "Pearson chi-square",
        Statistic = "Not estimated", df = "-", P_value = "-", Significance = "-",
        Note = "Chi-square failed because of sparse cells"
      )
    } else {
      low_expected <- any(tst$expected < 5, na.rm = TRUE)
      res[[length(res) + 1]] <- data.frame(
        Variable = v,
        Variable_Type = "Categorical / factor",
        Test = "Pearson chi-square",
        Statistic = paste0("Chi-square = ", round(as.numeric(tst$statistic), 3)),
        df = as.character(as.numeric(tst$parameter)),
        P_value = format_p_apa(tst$p.value),
        Significance = sig_label(tst$p.value),
        Note = ifelse(low_expected, "Some expected cells < 5; interpret cautiously", interpret_p(tst$p.value))
      )
    }
  }

  for (v in all_cont) {
    x <- safe_numeric(df[[v]])
    tmp <- data.frame(y = y_fac, x = x) %>% tidyr::drop_na()
    if (nrow(tmp) < 3 || dplyr::n_distinct(tmp$y) < 2 || dplyr::n_distinct(tmp$x) < 2) {
      res[[length(res) + 1]] <- data.frame(
        Variable = v, Variable_Type = "Numerical", Test = "One-way ANOVA",
        Statistic = "Not estimated", df = "-", P_value = "-", Significance = "-",
        Note = "Insufficient variation or sample size"
      )
      next
    }
    mod <- tryCatch(stats::aov(x ~ y, data = tmp), error = function(e) NULL)
    if (is.null(mod)) {
      res[[length(res) + 1]] <- data.frame(
        Variable = v, Variable_Type = "Numerical", Test = "One-way ANOVA",
        Statistic = "Not estimated", df = "-", P_value = "-", Significance = "-",
        Note = "ANOVA model failed"
      )
    } else {
      s <- summary(mod)[[1]]
      fval <- s$`F value`[1]
      pval <- s$`Pr(>F)`[1]
      df1 <- s$Df[1]
      df2 <- s$Df[2]
      res[[length(res) + 1]] <- data.frame(
        Variable = v,
        Variable_Type = "Numerical",
        Test = "One-way ANOVA",
        Statistic = paste0("F(", df1, ", ", df2, ") = ", round(fval, 3)),
        df = paste0(df1, ", ", df2),
        P_value = format_p_apa(pval),
        Significance = sig_label(pval),
        Note = interpret_p(pval)
      )
    }
  }

  if (length(res) == 0) {
    return(data.frame(Message = "Select categorical/factor or numerical predictors in the Data Setup tab."))
  }
  dplyr::bind_rows(res)
}

calc_inequality_apa_tables <- function(spatial_df, wealth_available = FALSE) {
  if (is.null(spatial_df)) return(list(summary = data.frame(Message = "Spatial data not available."), regional = data.frame()))
  df <- spatial_df %>% sf::st_drop_geometry() %>% dplyr::filter(!is.na(outcome))
  if (nrow(df) == 0) return(list(summary = data.frame(Message = "Outcome values not available."), regional = data.frame()))

  g_val <- ineq::ineq(df$outcome, type = "Gini")
  t_val <- ineq::ineq(df$outcome, type = "Theil")
  mu <- mean(df$outcome, na.rm = TRUE)
  sdv <- stats::sd(df$outcome, na.rm = TRUE)
  cv <- ifelse(mu == 0, NA, sdv / mu)
  min_v <- min(df$outcome, na.rm = TRUE)
  max_v <- max(df$outcome, na.rm = TRUE)
  ratio_v <- ifelse(min_v == 0, NA, max_v / min_v)
  range_v <- max_v - min_v

  conc <- NA_real_
  if (wealth_available && "bi_var" %in% names(df) && any(!is.na(df$bi_var))) {
    dfi <- df %>% dplyr::filter(!is.na(bi_var)) %>% dplyr::arrange(bi_var)
    if (nrow(dfi) > 1 && mean(dfi$outcome, na.rm = TRUE) != 0) {
      conc <- (2 / (nrow(dfi) * mean(dfi$outcome, na.rm = TRUE))) * sum(dfi$outcome * seq_len(nrow(dfi))) - 1 - (1 / nrow(dfi))
    }
  }

  summary_tbl <- data.frame(
    Metric = c("Gini coefficient", "Theil index", "Concentration index", "Mean regional outcome", "Standard deviation", "Coefficient of variation", "Minimum", "Maximum", "Absolute range", "Maximum/minimum ratio"),
    Estimate = c(g_val, t_val, conc, mu, sdv, cv, min_v, max_v, range_v, ratio_v),
    APA_Formatted = c(round(g_val, 4), round(t_val, 4), ifelse(is.na(conc), "-", round(conc, 4)), round(mu, 4), round(sdv, 4), ifelse(is.na(cv), "-", round(cv, 4)), round(min_v, 4), round(max_v, 4), round(range_v, 4), ifelse(is.na(ratio_v), "-", round(ratio_v, 4))),
    Interpretation = c(
      "Higher values indicate stronger geographic concentration.",
      "Higher values indicate stronger decomposable inequality.",
      "Positive values indicate concentration among higher-ranked areas on the selected ranking variable.",
      "Average outcome across mapped areas.",
      "Absolute dispersion across mapped areas.",
      "Relative dispersion around the mean.",
      "Lowest mapped regional estimate.",
      "Highest mapped regional estimate.",
      "Difference between highest and lowest mapped estimate.",
      "Relative gap between highest and lowest mapped estimate."
    )
  )

  regional_tbl <- df %>%
    dplyr::mutate(
      Region = Match_Key,
      Outcome = outcome,
      Outcome_Share_Percent = ifelse(sum(outcome, na.rm = TRUE) == 0, NA, outcome / sum(outcome, na.rm = TRUE) * 100),
      Rank_Low_to_High = rank(outcome, ties.method = "average"),
      Quintile = dplyr::ntile(outcome, 5),
      Gap_From_Mean = outcome - mu,
      Relative_Gap_Percent = ifelse(mu == 0, NA, (outcome - mu) / mu * 100)
    ) %>%
    dplyr::arrange(Outcome) %>%
    dplyr::mutate(Cumulative_Share_Percent = cumsum(Outcome_Share_Percent)) %>%
    dplyr::select(Region, n_obs, Outcome, Outcome_Share_Percent, Cumulative_Share_Percent, Rank_Low_to_High, Quintile, Gap_From_Mean, Relative_Gap_Percent)

  list(summary = summary_tbl, regional = regional_tbl)
}



# --- FINAL READY VERSION HELPERS: robust model tables and diagnostics ---
tidy_glmmTMB_fixed_apa <- function(model, exponentiate = TRUE) {
  if (is.null(model)) return(data.frame(Message = "Model was not run."))
  out <- tryCatch({
    sm <- summary(model)
    co <- sm$coefficients$cond
    if (is.null(co)) co <- sm$coefficients
    co <- as.data.frame(co)
    co$Term <- rownames(co)
    names(co) <- gsub(" ", "_", names(co), fixed = TRUE)
    est_col <- intersect(c("Estimate", "estimate"), names(co))[1]
    se_col  <- intersect(c("Std._Error", "Std. Error", "std.error"), names(co))[1]
    p_col   <- intersect(c("Pr(>|z|)", "Pr(>|t|)", "p.value"), names(co))[1]
    z_col   <- intersect(c("z_value", "z.value", "t_value", "t.value"), names(co))[1]
    est <- as.numeric(co[[est_col]])
    se <- if(!is.na(se_col)) as.numeric(co[[se_col]]) else rep(NA_real_, length(est))
    stat <- if(!is.na(z_col)) as.numeric(co[[z_col]]) else est / se
    pval <- if(!is.na(p_col)) as.numeric(co[[p_col]]) else 2 * (1 - pnorm(abs(stat)))
    lci <- est - 1.96 * se
    uci <- est + 1.96 * se
    if (exponentiate) {
      data.frame(
        Term = co$Term,
        Estimate = round(est, 4),
        AOR = round(exp(est), 4),
        LCI = round(exp(lci), 4),
        UCI = round(exp(uci), 4),
        Statistic = round(stat, 3),
        P_Value = sapply(pval, format_p_apa),
        Significance = sapply(pval, sig_label),
        check.names = FALSE
      )
    } else {
      data.frame(
        Term = co$Term,
        Estimate = round(est, 4),
        LCI = round(lci, 4),
        UCI = round(uci, 4),
        Statistic = round(stat, 3),
        P_Value = sapply(pval, format_p_apa),
        Significance = sapply(pval, sig_label),
        check.names = FALSE
      )
    }
  }, error = function(e) {
    data.frame(Message = paste("Could not tidy model:", e$message))
  })
  out
}

safe_spatial_weights <- function(df, nb_type = "Queen's Contiguity") {
  nb <- tryCatch({
    if(nb_type == "Queen's Contiguity") {
      spdep::poly2nb(df, queen = TRUE)
    } else if(nb_type == "K-Nearest Neighbors (K=3)") {
      spdep::knn2nb(spdep::knearneigh(sf::st_coordinates(sf::st_centroid(df)), k = 3))
    } else {
      spdep::knn2nb(spdep::knearneigh(sf::st_coordinates(sf::st_centroid(df)), k = 5))
    }
  }, error = function(e) spdep::knn2nb(spdep::knearneigh(sf::st_coordinates(sf::st_centroid(df)), k = 1)))
  if(any(spdep::card(nb) == 0)) {
    nb <- spdep::knn2nb(spdep::knearneigh(sf::st_coordinates(sf::st_centroid(df)), k = 1))
  }
  lw <- spdep::nb2listw(nb, style = "W", zero.policy = TRUE)
  list(nb = nb, lw = lw)
}

ui <- dashboardPage(
  dark = NULL,
  header = dashboardHeader(title = "HawaSpatial Pro"),
  sidebar = dashboardSidebar(
    sidebarMenu(
      menuItem("1. Introduction", tabName = "intro", icon = icon("info-circle")),
      menuItem("2. Data Setup & Equalizer", tabName = "setup", icon = icon("database")),
      menuItem("3. 17 Spatial Methods", tabName = "esda", icon = icon("map-marked-alt")),
      menuItem("4. Bivariate Spatial", tabName = "bivariate", icon = icon("layer-group")),
      menuItem("5. Advanced Spatial", tabName = "advanced", icon = icon("magic")),
      menuItem("6. Multilevel Models", tabName = "me_ml", icon = icon("brain")),
      menuItem("7. Spatial Random Effects", tabName = "me_spatial", icon = icon("map-marker-alt")),
      menuItem("8. Small Area Estimation", tabName = "sae_tab", icon = icon("chart-area")),
      menuItem("9. Spatial Econometrics", tabName = "spatial_econ", icon = icon("calculator")),
      menuItem("10. Spatial Inequality", tabName = "inequality", icon = icon("balance-scale")),
      menuItem("11. Regional Decomposition", tabName = "spatial_decomp", icon = icon("scissors")),
      menuItem("12. Reporting", tabName = "reporting", icon = icon("file-export")),
      menuItem("13. Spatiotemporal Dynamics", tabName = "st_dynamics", icon = icon("history")),
      menuItem("14. Spatiotemporal Decomp", tabName = "decomp", icon = icon("scissors")),
      menuItem("15. Joint Spatial Modeling", tabName = "joint_mod", icon = icon("project-diagram")),
      menuItem("16. Joint Regional Decomp", tabName = "joint_decomp", icon = icon("scissors")),
      menuItem("17. Joint Spatiotemporal Decomp", tabName = "joint_st_decomp", icon = icon("scissors")),
      menuItem("18. Paper 11: Service Deserts", tabName = "paper11_access", icon = icon("hospital")),
      menuItem("19. Paper 12: Local Determinants", tabName = "paper12_gwr", icon = icon("map-location-dot")),
      menuItem("20. Paper 13: Priority Index", tabName = "paper13_priority", icon = icon("bullseye")),
      menuItem("21. Paper 14: Cluster Detection", tabName = "paper14_scan", icon = icon("search-location")),
      menuItem("22. Paper 15: Spatial Equity", tabName = "paper15_equity", icon = icon("scale-balanced"))
    )
  ),
  body = dashboardBody(
    tags$head(tags$style(HTML(custom_theme))),
    tabItems(
      tabItem(tabName = "intro",
              jumbotron(
                title = "HAWA Spatial Intelligence",
                lead = "Holistic & Areal Weighted Analysis (HAWA) for Global Development",
                status = "success",
                "HawaSpatial Pro is a professional-grade, multi-sectoral platform designed to bridge the gap
                between complex spatial modeling and actionable policy insights. It integrates
                advanced computational statistics with Geographic Information Systems (GIS) to provide
                a robust framework for sub-national monitoring of development indicators. The platform
                supports comprehensive analysis across seven domains: Exploratory Spatial Data Analysis,
                Bivariate Spatial Association, Spatiotemporal Dynamics, Small Area Estimation (SAE),
                Spatial Econometrics, Spatial Inequality Metrics, and Hierarchical Multilevel Modeling.
                The platform is optimized for large-scale, complex household surveys including DHS, MICS,
                IHBS, MIS, SPA, and LSMS, supporting both single-survey cross-sectional analysis and
                multiple survey waves for longitudinal tracking."
              ),
              fluidRow(
                box(title = "Core Analytical Framework", width = 12, status = "primary", solidHeader = TRUE,
                    column(3,
                           h4(icon("chart-pie"), "1. Descriptive"),
                           p("Establish baselines and measure concentration:"),
                           tags$ul(
                             tags$li(strong("Data Integration:"), " Features a 'Spatial Equalizer' for seamless merging of complex survey data with administrative shapefiles."),
                             tags$li(strong("Weighted Stats:"), " Generates weighted frequency, %, and Mean (SD) characteristics of respondents."),
                             tags$li(strong("Spatial Inequality:"), " Calculates Gini Coefficient and Theil Index with Lorenz Curve visualization to measure geographic concentration.")
                           )),
                    column(3,
                           h4(icon("microscope"), "2. Diagnostic"),
                           p("Identify patterns and drivers of change:"),
                           tags$ul(
                             tags$li(strong("17 Spatial Methods:"), " Includes Global/Local Moran's I, LISA Cluster Maps, Getis-Ord Gi*, Spatial Correlograms, and Join Count Analysis."),
                             tags$li(strong("Bivariate & Advanced:"), " Measures spatial covariance via Bivariate Moran's I and Ordinary Kriging surface interpolation."),
                             tags$li(strong("Dynamics:"), " Analyzes trends via Rate of Change maps, LISA Transitions, and Oaxaca-Blinder style Multivariate Decomposition.")
                           )),
                    column(3,
                           h4(icon("brain"), "3. Predictive"),
                           p("Stabilize estimates and model risk:"),
                           tags$ul(
                             tags$li(strong("Multilevel Modeling:"), " Runs 4-Stage Hierarchical Models (Null, Individual, Community, Full) via the glmmTMB engine."),
                             tags$li(strong("Small Area Estimation:"), " Implements Fay-Herriot area-level models to produce stabilized sub-national EBLUPs."),
                             tags$li(strong("Spatial Econometrics:"), " Features Spatial Autoregressive (SAR) and Spatial Error Models (SEM) to quantify spillover effects.")
                           )),
                    column(3,
                           h4(icon("bullseye"), "4. Prescriptive"),
                           p("Evidence-based policy and reporting:"),
                           tags$ul(
                             tags$li(strong("Target Gap Analysis:"), " Identifies 'SDG Target Gaps' to enable precision targeting of hotspots."),
                             tags$li(strong("Random Effects:"), " Maps unexplained variation to identify unobserved geographic risk factors."),
                             tags$li(strong("Automated Reporting:"), " Generates formal publication-ready tables and reports in Word, PDF, and HTML formats.")
                           ))
                )
              ),

              # --- NEW SECTION: THE 17-TAB INTEGRATED WORKFLOW ---
              fluidRow(
                box(title = "The 22-Tab Integrated Workflow", width = 12, status = "info", solidHeader = TRUE,
                    p("The platform is organized into 22 specialized tabs that guide the user through a complete scientific discovery process:"),
                    column(4,
                           tags$ol(
                             tags$li(strong("Introduction:"), " Project overview and methodological framework."),
                             tags$li(strong("Data Setup:"), " Smart-matching of survey data and shapefiles."),
                             tags$li(strong("17 Spatial Methods:"), " Comprehensive ESDA diagnostic suite."),
                             tags$li(strong("Bivariate Spatial:"), " Spatial covariance and Bivariate Moran's I."),
                             tags$li(strong("Advanced Spatial:"), " Ordinary Kriging surface interpolation."),
                             tags$li(strong("Multilevel Models:"), " 4-stage hierarchical linear modeling.")
                           )),
                    column(4,
                           tags$ol(start = 7,
                                   tags$li(strong("Spatial Random Effects:"), " Mapping unexplained cluster variation."),
                                   tags$li(strong("Small Area Estimation:"), " Fay-Herriot stabilized EBLUPs."),
                                   tags$li(strong("Spatial Econometrics:"), " Multivariate SAR and SEM models."),
                                   tags$li(strong("Spatial Inequality:"), " Gini, Theil, and Concentration indices."),
                                   tags$li(strong("Regional Decomposition:"), " Oaxaca-Blinder geographic gap analysis."),
                                   tags$li(strong("Reporting:"), " Publication-ready automated exports.")
                           )),
                    column(4,
                           tags$ol(start = 13,
                                   tags$li(strong("Spatiotemporal Dynamics:"), " Trend analysis and LISA transitions."),
                                   tags$li(strong("Spatiotemporal Decomp:"), " Oaxaca-Blinder drivers of change over time."),
                                   tags$li(strong("Joint Spatial Modeling:"), " Shared determinants of syndemic burdens."),
                                   tags$li(strong("Joint Regional Decomp:"), " Explaining geographic gaps in joint outcomes."),
                                   tags$li(strong("Joint Spatiotemporal Decomp:"), " Evolution of the syndemic burden over time."),
                                   tags$li(strong("Service Deserts:"), " Distance/accessibility gaps and underserved areas."),
                                   tags$li(strong("Local Determinants:"), " Spatially varying predictor effects."),
                                   tags$li(strong("Priority Index:"), " Composite spatial targeting score."),
                                   tags$li(strong("Cluster Detection:"), " High-burden scan-style cluster summaries."),
                                   tags$li(strong("Spatial Equity:"), " Equity gap analysis by vulnerable group.")
                           ))
                )
              ),

              # --- NEW SECTION: THE 9-PAPER PUBLICATION SUITE ---
              fluidRow(
                box(title = "The 15-Paper Publication Suite", width = 12, status = "danger", solidHeader = TRUE,
                    p("HawaSpatial Pro is designed to generate the core results for fifteen distinct, high-impact scientific papers from a single data run:"),
                    column(4,
                           tags$ul(
                             tags$li(strong("Paper 1:"), " Small Area Estimation (SAE) for precision sub-national monitoring."),
                             tags$li(strong("Paper 2:"), " Multivariate Spatial Econometrics and Spillover Analysis."),
                             tags$li(strong("Paper 3:"), " Spatial Inequality and SDG 10 Monitoring.")
                           )),
                    column(4,
                           tags$ul(
                             tags$li(strong("Paper 4:"), " Multilevel Determinants and Variance Decomposition."),
                             tags$li(strong("Paper 5:"), " Decomposing the Drivers of Regional Inequality."),
                             tags$li(strong("Paper 6:"), " Spatiotemporal Dynamics and Drivers of Change.")
                           )),
                    column(4,
                           tags$ul(
                             tags$li(strong("Paper 7:"), " Spatiotemporal Decomposition of Change."),
                             tags$li(strong("Paper 8:"), " Joint Spatial Modeling and Shared Determinants (Syndemics)."),
                             tags$li(strong("Paper 9:"), " Decomposing the Regional Gap in the Syndemic Burden."),
                             tags$li(strong("Paper 10:"), " Spatiotemporal Drivers of the Syndemic Evolution."),
                             tags$li(strong("Paper 11:"), " Spatial Accessibility and Service Desert Analysis."),
                             tags$li(strong("Paper 12:"), " Geographically Weighted / Local Determinants Analysis."),
                             tags$li(strong("Paper 13:"), " Spatial Priority Index and Targeting Analysis."),
                             tags$li(strong("Paper 14:"), " Spatial Scan and Cluster Detection Analysis."),
                             tags$li(strong("Paper 15:"), " Spatial Equity and Environmental Justice Analysis.")
                           ))
                )
              ),

              fluidRow(
                box(title = "Universal SDG Alignment", width = 12, status = "success", solidHeader = TRUE,
                    p("HawaSpatial Pro is engineered to support the monitoring and achievement of all 17 United Nations Sustainable Development Goals (SDGs).
                      By identifying geographic 'Hot Spots' and 'Inequality Gaps', the platform enables precision targeting for:"),
                    column(4,
                           tags$ul(
                             tags$li(strong("SDG 1-5:"), " Poverty, Hunger, Health, Education, & Gender Equality"),
                             tags$li(strong("SDG 10:"), " Reduced Inequalities (Geographic Gini & Theil Indices)")
                           )),
                    column(4,
                           tags$ul(
                             tags$li(strong("SDG 6-9:"), " Water, Energy, Decent Work, & Infrastructure"),
                             tags$li(strong("SDG 11-15:"), " Sustainable Cities, Climate Action, & Life on Land")
                           )),
                    column(4,
                           tags$ul(
                             tags$li(strong("SDG 16-17:"), " Peace, Justice, & Global Partnerships for the Goals"),
                             tags$li(strong("Outcome Flexibility:"), " Supports Binary, Count, Continuous, Ordinal, Multinomial, and Fractional outcomes.")
                           ))
                )
              ),
              fluidRow(
                box(title = "Project Overview & Data Sources", width = 6, status = "primary", solidHeader = TRUE,
                    p("The platform is optimized for large-scale, complex household surveys including:"),
                    tags$ul(
                      tags$li(strong("DHS:"), " Demographic and Health Surveys (USAID)"),
                      tags$li(strong("MICS:"), " Multiple Indicator Cluster Surveys (UNICEF)"),
                      tags$li(strong("IHBS/LSMS:"), " Integrated Household Budget & Living Standards Measurement Surveys"),
                      tags$li(strong("SPA/MIS:"), " Service Provision Assessment & Malaria Indicator Surveys")
                    ),
                    p("It features a built-in 'Spatial Equalizer' to resolve naming inconsistencies between survey datasets and administrative shapefiles (Admin1/Admin2).")
                ),
                box(title = "Technical Engine & R Packages", width = 6, status = "info", solidHeader = TRUE,
                    p("HawaSpatial Pro leverages a sophisticated ecosystem of R libraries:"),
                    tags$ul(
                      tags$li(strong("glmmTMB, performance & survey:"), " used for Hierarchical Models, ICC/PCV calculation, and Design-Based Inference."),
                      tags$li(strong("sae & spatialreg:"), " used for Small Area Estimation (Fay-Herriot) and Spatial Econometrics (SAR/SEM/Impacts)."),
                      tags$li(strong("ineq & biscale:"), " used for Gini, Theil, and Bivariate mapping."),
                      tags$li(strong("sf, spdep & gstat:"), " Simple Features, Spatial Dependence, and Geostatistics for weights and Kriging.")
                    )
                )
              )
      ),
      tabItem(tabName = "setup",
              fluidRow(
                box(title = "1. Data Ingestion", width = 4, status = "primary",
                    fileInput("data_file", "Upload Survey Data (.dta, .sav, .csv)"),
                    fileInput("shp_file", "Upload Shapefile Components", multiple = TRUE),
                    selectInput("outcome_type", "Outcome Type:",
                                choices = c("Binary", "Count", "Continuous", "Ordinal", "Multinomial", "Fractional")),
                    checkboxInput("is_dhs_wgt", "Apply DHS Weight (v005/1,000,000)", TRUE),
                    actionButton("process_data", "Initialize & Equalize Names", class = "btn-success btn-block")
                ),
                box(title = "2. Core Mapping Variables", width = 8, status = "primary",
                    fluidRow(
                      column(6, selectInput("v_dv", "Outcome (dv):", choices = NULL),
                             selectInput("v_reg", "Region Variable:", choices = NULL)),
                      column(6, selectInput("v_cluster", "Cluster ID (V001):", choices = NULL),
                             selectInput("v_wgt", "Weight (V005):", choices = NULL)),
                      selectInput("v_strata", "Strata Variable (V023):", choices = NULL),
                    ),
                    hr(),
                    fluidRow(
                      column(6, selectInput("v_indiv_cont", "Individual Continuous (L1):", choices = NULL, multiple = TRUE),
                             selectInput("v_indiv_fact", "Individual Factors (L1):", choices = NULL, multiple = TRUE)),
                      column(6, selectInput("v_comm_cont", "Community Continuous (L2):", choices = NULL, multiple = TRUE),
                             selectInput("v_comm_fact", "Community Factors (L2):", choices = NULL, multiple = TRUE))
                    ),
                    hr(),
                    # --- NEW: TIME VARIABLE INPUT (Optional) ---
                    selectInput("v_time", "Time/Wave Variable (for Spatiotemporal):", choices = NULL),
                    selectInput("v_wealth", "Bivariate X Variable (e.g. Wealth):", choices = NULL),
                    selectInput("v_risk_group", "Variable for Regional Groups (e.g. Risk Level):", choices = NULL),
                    uiOutput("sparsity_alert"))
              )
      ),
      tabItem(tabName = "esda",
              fluidRow(
                box(title = "ESDA Diagnostic Selector", width = 3, status = "success",
                    selectInput("esda_method", "Select Method:",
                                choices = c("1. Connectivity Graph", "2. Basic Choropleth", "3. EBS Smoothing",
                                            "4. Quantile Map", "5. Global Moran's I", "6. Moran Scatterplot",
                                            "7. LISA Cluster Map", "8. Getis-Ord Gi*", "9. Join Count Analysis",
                                            "10. Hinge Map", "11. Spatial Correlogram", "12. Regime Density",
                                            "13. Spatial Lag Map", "14. SDG Target Gap", "15. Local Geary's C",
                                            "16. Conditional Plot", "17. Global Geary's C"))
                ),
                box(title = "Spatial Diagnostic Output", width = 9, plotOutput("diag_plot", height = "600px"))
              ),
              fluidRow(
                box(title = "Table 1: Characteristics of Respondents (Weighted)", width = 12, status = "primary", solidHeader = TRUE,
                    helpText("Displays Weighted Frequency and % for Categorical variables, and Weighted Mean (SD) for Continuous variables."),
                    tableOutput("desc_stats_table")
                )
              )
              # --------------------------------------
      ),

      tabItem(tabName = "bivariate",
              fluidRow(
                box(title = "Bivariate Controls", width = 4, status = "warning",
                    actionButton("run_bi", "Generate Bivariate Map", class = "btn-success"),
                    hr(),
                    h4("Bivariate Global Moran's I Test"),
                    verbatimTextOutput("bi_moran_test") # New output for Bivariate Moran
                ),
                box(title = "Bivariate Spatial Relationship", width = 8, plotOutput("bi_plot", height = "600px"))
              ),

              fluidRow(
                box(title = "Multicollinearity Check: Covariate Correlation Matrix", width = 12, status = "primary", solidHeader = TRUE,
                    plotOutput("cov_corr_plot", height = "600px"),
                    helpText("Visualizes correlations between all selected Individual and Community level independent variables. Darker colors indicate stronger correlation."))
              )
              # ------------------------
      ),

      tabItem(tabName = "advanced",
              fluidRow(
                box(title = "Ordinary Kriging (Surface Interpolation)", width = 6, plotOutput("krig_plot")),
                box(title = "Spatial Lag Model Residuals", width = 6, plotOutput("resid_plot"))
              ),
              fluidRow(
                box(title = "APA Table: Bivariate Screening Results", width = 12, status = "success", solidHeader = TRUE,
                    helpText("Categorical/factor predictors are tested using Pearson chi-square. Numerical predictors are tested using one-way ANOVA across outcome groups."),
                    tableOutput("advanced_biv_table"))
              )
      ),

      tabItem(tabName = "me_ml",
              fluidRow(
                box(width = 3, title = "Multilevel Engine", status = "info",
                    actionButton("run_freq", "Run 4-Stage Hierarchical Model", class = "btn-info btn-block"),
                    hr(),
                    helpText("Generates comparative models M0-M3 for publication.")
                ),
                box(title = "Full Model: Fixed Effect Estimates", width = 9, plotOutput("freq_forest_plot"))
              ),
              fluidRow(
                box(title = "Formal Publication Table: Side-by-Side Model Comparison", width = 12,
                    status = "success", solidHeader = TRUE,
                    helpText("Includes Models 0, 1, 2, and 3. Estimates are reported as Adjusted Odds Ratios (AOR) for Binary/Count outcomes."),
                    tableOutput("ml_comp_table"))
              ),
              fluidRow(
                box(title = "Model Fitness & Variance Metrics", width = 12, plotOutput("freq_metrics_plot"))
              )
      ),
      tabItem(tabName = "me_spatial",
              fluidRow(
                box(title = "Unexplained Spatial Variation (Random Effects)", width = 12, status = "info",
                    plotOutput("freq_map", height = "700px"))
              )
      ),
      tabItem(tabName = "sae_tab",
              fluidRow(
                box(title = "Small Area Estimation (Fay-Herriot Model)", width = 4, status = "primary",
                    actionButton("run_sae", "Run SAE Model", class = "btn-success btn-block"),
                    hr(),
                    helpText("This module uses the Fay-Herriot area-level model to stabilize regional estimates. It combines direct survey estimates with a regression model (using your Bivariate X variable) to reduce variance in regions with small sample sizes.")
                ),
                box(title = "SAE Stabilized Estimates (EBLUPs)", width = 8, plotOutput("sae_map", height = "600px"))
              ),
              fluidRow(
                box(title = "SAE Model Comparison: Direct vs. EBLUP", width = 12,
                    status = "info", solidHeader = TRUE,
                    tableOutput("sae_table"))
              ),
              fluidRow(
                box(title = "SAE Model Diagnostics", width = 4, status = "info",
                    tableOutput("sae_diag_table")),
                box(title = "Shrinkage Plot (Direct vs. EBLUP)", width = 8,
                    plotOutput("sae_shrinkage_plot"))
              )
      ),
      tabItem(tabName = "inequality",
              fluidRow(
                box(title = "SDG 10: Spatial Inequality Metrics", width = 4, status = "danger",
                    infoBoxOutput("gini_box", width = 12),
                    infoBoxOutput("theil_box", width = 12),
                    infoBoxOutput("conc_box", width = 12),
                    hr(),
                    helpText("The Gini Coefficient (0 to 1) measures geographic concentration. 1 indicates total inequality. The Theil Index decomposes inequality and is sensitive to outliers.")
                ),
                box(title = "Lorenz Curve (Geographic Distribution)", width = 8,
                    status = "primary", solidHeader = TRUE,
                    plotOutput("lorenz_plot", height = "500px"),
                    helpText("The Lorenz Curve visualizes the cumulative distribution of the outcome across regions. The further the curve is from the diagonal line, the higher the geographic inequality."))
              ),
              fluidRow(
                box(title = "APA Table: Spatial Inequality Summary", width = 12, status = "success", solidHeader = TRUE,
                    tableOutput("ineq_apa_table"))
              ),
              fluidRow(
                box(title = "Detailed Regional Inequality Ranking", width = 12, status = "primary", solidHeader = TRUE,
                    helpText("This table ranks mapped areas by outcome level, share of total burden, cumulative share, quintile, and gap from the mean."),
                    tableOutput("ineq_region_detail_table"))
              )
      ),
      tabItem(tabName = "spatial_econ",
              fluidRow(
                box(title = "1. Model Selection (Lagrange Multiplier Tests)", width = 12, status = "warning",
                    verbatimTextOutput("econ_lm_tests"),
                    helpText("LM tests help choose between SAR and SEM. If LM-Lag is significant, use SAR. If LM-Error is significant, use SEM."))
              ),
              fluidRow(
                box(title = "2. Multivariate Spatial Regression", width = 4, status = "primary",
                    selectInput("econ_nb_type", "Spatial Weights Matrix:",
                                choices = c("Queen's Contiguity", "K-Nearest Neighbors (K=3)", "K-Nearest Neighbors (K=5)")),
                    selectInput("econ_model_type", "Select Model Type:",
                                choices = c("Spatial Lag Model (SAR)", "Spatial Error Model (SEM)")),
                    actionButton("run_econ", "Run Multivariate Regression", class = "btn-success btn-block"),
                    hr(),
                    helpText(strong("Note on Interpretation:"), "This is an aggregate-level model. Results for variables like 'Education' may differ from individual-level models due to the Ecological Paradox (e.g., individual protection vs. regional backlash).")
                ),
                box(title = "Model Coefficients", width = 8, status = "info", tableOutput("econ_coef_table"))
              ),
              fluidRow(
                box(title = "3. Spillover Impacts (Direct, Indirect, Total)", width = 12, status = "success",
                    tableOutput("econ_impacts_table"))
              ),
              fluidRow(
                box(title = "4. APA Model Fit and Spatial Diagnostic Table", width = 6, status = "info", solidHeader = TRUE,
                    tableOutput("econ_fit_table")),
                box(title = "5. Spatial Weights Summary", width = 6, status = "info", solidHeader = TRUE,
                    tableOutput("econ_weights_table"))
              ),
              fluidRow(
                box(title = "6. Regional Observed, Fitted, and Residual Values", width = 12, status = "primary", solidHeader = TRUE,
                    tableOutput("econ_resid_table"))
              )
      ),
      tabItem(tabName = "reporting",
              fluidRow(
                box(title = "Export Publication Report", width = 6, status = "success", solidHeader = TRUE,
                    textInput("rep_title", "Report Title:", "HAWA Spatial Intelligence Publication Report"),
                    radioButtons("format", "Format:", choices = c("Word" = "docx", "PDF" = "pdf", "HTML" = "html")),
                    downloadButton("dl_report", "Download Full Report")
                ),
                box(title = "How to Cite HawaSpatial Pro", width = 6, status = "info", solidHeader = TRUE,
                    p("If you use HawaSpatial Pro in your research, please cite it as:"),
                    tags$blockquote(
                      "Muse, A. H. (2026). HawaSpatial: Holistic and Areal Weighted Analysis for Global Development (Version 0.1.10). Zenodo. https://doi.org/10.5281/zenodo.18290101"
                    ),
                    hr(),
                    h5("BibTeX for LaTeX users:"),
                    tags$pre(
                      "@manual{muse2026hawa,\n  title = {HawaSpatial: Holistic and Areal Weighted Analysis for Global Development},\n  author = {Muse, Abdisalam Hassan},\n  year = {2026},\n  note = {R package version 0.1.10},\n  doi = {10.5281/zenodo.18290101},\n  url = {https://doi.org/10.5281/zenodo.18290101}\n}"
                    )
                )
              )
      ),
      # --- NEW: SPATIOTEMPORAL DYNAMICS TAB ---
      tabItem(tabName = "decomp",
              fluidRow(
                box(title = "Spatiotemporal Decomposition Controls", width = 4, status = "primary",
                    helpText("Compare two time periods (Waves) to see what drove the change."),
                    uiOutput("wave_selector_1"), uiOutput("wave_selector_2"),
                    actionButton("run_decomp", "Run Spatiotemporal Decomposition", class = "btn-success btn-block")
                ),
                box(title = "Table 4: Summary Decomposition (Time Gap Partitioning)", width = 8, status = "info",
                    tableOutput("st_decomp_summary_table")),
                box(title = "Table 2: Descriptive Profile of Waves (Weighted %)", width = 12, status = "primary",
                    tableOutput("st_table2"))
              ),
              fluidRow(
                box(title = "Table 3: Outcome Rates by Characteristics (Wave Comparison)", width = 12,
                    status = "primary", tableOutput("st_table3"))
              ),
              fluidRow(
                box(title = "Figure 5: Visualizing Drivers of Change (Forest Plot)", width = 7,
                    plotOutput("st_figure5_plot", height = "500px")),
                box(title = "Summary: Percentage Contribution to Change", width = 5,
                    plotOutput("st_pct_plot")),
                box(title = "Table 5/6: Detailed Variable Contributions (Time)", width = 5,
                    tableOutput("st_decomp_detail_table"))
              )
      ),
      tabItem(tabName = "st_dynamics",
              fluidRow(
                box(title = "4.1 Global Spatial Trend Over Time", width = 12, status = "primary",
                    uiOutput("st_moran_trend_ui"),
                    hr(),
                    h4("APA Table: Moran's I Trend and P-values"),
                    tableOutput("st_moran_trend_table")) # Use UI output for conditional rendering
              ),
              fluidRow(
                box(title = "4.2 Outcome Rate by Time Period (Small Multiples)", width = 12, status = "info",
                    plotOutput("st_small_multiples", height = "800px"))
              ),
              fluidRow(
                box(title = "4.3 Outcome Rate of Change (Min to Max Time)", width = 6, status = "info",
                    uiOutput("st_rate_change_ui")),
                box(title = "4.4 LISA Cluster Transition Map", width = 6, status = "info",
                    uiOutput("st_lisa_transition_ui"))
              ),
              fluidRow(
                box(title = "APA Table: Rate of Change by Region", width = 6, status = "primary", solidHeader = TRUE,
                    tableOutput("st_rate_change_table")),
                box(title = "APA Table: LISA Transition Classification", width = 6, status = "primary", solidHeader = TRUE,
                    tableOutput("st_lisa_transition_table"))
              ),
              hr(),
              h3(icon("Microscope"), "4.5 Predictive Spatio-Temporal Interaction Model"),
              p("This model (Bernardinelli-style) partitions risk into a global trend and regional deviations."),
              fluidRow(
                box(title = "Interaction Model Controls", width = 4, status = "danger", solidHeader = TRUE,
                    helpText("Formula: dv ~ predictors + time + (1|Region) + (1|Year) + (1|Region:Year)"),
                    actionButton("run_st_int", "Run Interaction Model", class = "btn-danger btn-block"),
                    hr(),
                    tableOutput("st_int_metrics")
                ),
                box(title = "Map: Emerging Hotspots (Interaction Effects)", width = 8, status = "danger",
                    plotOutput("st_int_map", height = "500px"),
                    helpText("Red = Emerging Hotspots (Risk rising faster than national average). Blue = Improving faster than average."))
              ),
              fluidRow(
                box(title = "Model Coefficients (AOR)", width = 12, status = "info",
                    tableOutput("st_int_coef_table"))
              ),
              fluidRow(
                box(title = "APA Table: Emerging Hotspot Ranking", width = 12, status = "danger", solidHeader = TRUE,
                    tableOutput("st_hotspot_table"))
              )
      ),
      tabItem(tabName = "spatial_decomp",
              fluidRow(
                box(title = "Regional Decomposition Controls", width = 4, status = "danger",
                    checkboxInput("auto_lisa", "Auto-select Hot/Cold Spots (via LISA)", FALSE),
                    uiOutput("risk_selector_1"), uiOutput("risk_selector_2"),
                    actionButton("run_spatial_decomp", "Run Regional Decomposition", class = "btn-danger btn-block")
                ),
                box(title = "Table 4: Summary Decomposition (Gap Partitioning)", width = 8, status = "info",
                    tableOutput("spatial_decomp_summary_table")),
                box(title = "Table 2: Descriptive Profile of Regions (Weighted %)", width = 12, status = "primary",
                    tableOutput("spatial_table2"))
              ),
              fluidRow(
                box(title = "Table 3: Outcome Rates by Characteristics (Group Comparison)", width = 12,
                    status = "primary", tableOutput("spatial_table3"))
              ),
              fluidRow(
                box(title = "Figure 5: Visualizing Drivers of Inequality (Forest Plot)", width = 7,
                    plotOutput("spatial_figure5_plot", height = "500px")),
                box(title = "Summary: Percentage Contribution to Inequality", width = 5, plotOutput("spatial_pct_plot")),
                box(title = "Table 5/6: Detailed Variable Contributions", width = 5,
                    tableOutput("spatial_decomp_detail_table"))
              )
      ),
      tabItem(tabName = "joint_mod",
              fluidRow(
                box(title = "Joint Modeling Controls", width = 4, status = "primary",
                    helpText("Model two outcomes simultaneously to find shared determinants and spatial clustering."),
                    selectInput("joint_dv1", "Outcome 1 (e.g. Malnutrition):", choices = NULL),
                    selectInput("joint_dv2", "Outcome 2 (e.g. Food Insecurity):", choices = NULL),
                    actionButton("run_joint", "Run Joint Analysis", class = "btn-success btn-block")
                ),
                box(title = "1. Shared Variance & Global Moran's I", width = 8, status = "info",
                    fluidRow(
                      column(6, infoBoxOutput("joint_corr_box", width = 12)),
                      column(6, infoBoxOutput("joint_moran_box", width = 12))
                    ),
                    verbatimTextOutput("joint_moran_text")
                )
              ),
              fluidRow(
                box(title = "2. Joint Forest Plot: Shared Determinants", width = 7,
                    plotOutput("joint_forest_plot", height = "500px"),
                    helpText("Points show the effect of predictors on both outcomes. Overlapping points suggest shared drivers.")),
                box(title = "3. LISA Cluster Map: Joint Hotspots", width = 5,
                    plotOutput("joint_lisa_map", height = "500px"),
                    helpText("Red areas (High-High) are regions where BOTH conditions significantly cluster together."))
              ),
              fluidRow(
                box(title = "4. Joint Model Coefficients Table (AOR)", width = 12,
                    tableOutput("joint_coef_table"))
              )
      ),
      tabItem(tabName = "joint_decomp",
              fluidRow(
                box(title = "Joint Decomposition Controls", width = 4, status = "danger", solidHeader = TRUE,
                    helpText("Explain the GAP in the JOINT burden between two regional groups."),
                    selectInput("jd_dv1", "Outcome 1 (e.g. Stunting):", choices = NULL),
                    selectInput("jd_dv2", "Outcome 2 (e.g. Wasting):", choices = NULL),
                    selectInput("jd_group", "Regional Grouping Variable:", choices = NULL),
                    uiOutput("jd_risk_selector_1"),
                    uiOutput("jd_risk_selector_2"),
                    actionButton("run_joint_decomp", "Run Joint Decomposition", class = "btn-danger btn-block"),
                    hr(),
                    uiOutput("jd_status_ui")
                ),
                box(title = "Table 4: Summary of Joint Gap Partitioning", width = 8, status = "info", solidHeader = TRUE,
                    helpText("Partitioning the total regional gap into Explained (Endowment) and Unexplained (Coefficient) components."),
                    tableOutput("jd_summary_table"))
              ),

              fluidRow(
                box(title = "Table 2: Descriptive Profile of Joint Groups (Weighted %)", width = 12, status = "primary",
                    helpText("Comparison of predictor distributions between the two selected regional groups."),
                    tableOutput("jd_table2"))
              ),

              fluidRow(
                box(title = "Table 3: Joint Outcome Rates by Characteristics", width = 12, status = "primary",
                    helpText("The prevalence of the joint burden across different categories within each region."),
                    tableOutput("jd_table3"))
              ),

              fluidRow(
                box(title = "Figure 5: Drivers of the Joint Burden Gap (Forest Plot)", width = 7, status = "success",
                    helpText("Visualizing which variables contribute most to the explained (Endowment) portion of the gap."),
                    plotOutput("jd_forest_plot", height = "500px")),
                box(title = "Summary: Percentage Contribution to Joint Gap", width = 5, status = "success",
                    plotOutput("jd_pct_plot", height = "500px"))
              ),

              fluidRow(
                box(title = "Table 5: Detailed Aggregated Variable Contributions", width = 12, status = "info",
                    helpText("Numerical breakdown of how each aggregated variable contributes to the Endowment and Coefficient effects."),
                    tableOutput("jd_detail_table"))
              )
      ),
      tabItem(tabName = "joint_st_decomp",
              fluidRow(
                box(title = "Joint Spatiotemporal Controls", width = 4, status = "success", solidHeader = TRUE,
                    helpText("Explain the CHANGE in the JOINT burden between two time periods (Waves)."),
                    selectInput("jsd_dv1", "Outcome 1:", choices = NULL),
                    selectInput("jsd_dv2", "Outcome 2:", choices = NULL),
                    uiOutput("jsd_wave_selector_1"),
                    uiOutput("jsd_wave_selector_2"),
                    actionButton("run_jsd", "Run Spatiotemporal Joint Decomp", class = "btn-success btn-block"),
                    hr(),
                    uiOutput("jsd_status_ui")
                ),
                box(title = "Table 4: Summary of Joint Change Partitioning", width = 8, status = "info",
                    tableOutput("jsd_summary_table"))
              ),
              fluidRow(
                box(title = "Table 2: Descriptive Profile of Waves (Weighted %)", width = 12, status = "primary",
                    tableOutput("jsd_table2"))
              ),
              fluidRow(
                box(title = "Figure 5: Drivers of Joint Change (Forest Plot)", width = 7, status = "success",
                    plotOutput("jsd_forest_plot", height = "500px")),
                box(title = "Contribution Share (%)", width = 5, status = "success",
                    plotOutput("jsd_pct_plot", height = "500px"))
              ),
              fluidRow(
                box(title = "Table 5: Detailed Aggregated Contributions (Time)", width = 12, status = "info",
                    tableOutput("jsd_detail_table"))
              )
      ),
      tabItem(tabName = "paper11_access",
              fluidRow(
                box(title = "Paper 11 Controls: Spatial Accessibility and Service Deserts", width = 4, status = "primary", solidHeader = TRUE,
                    helpText("Optional: upload a CSV of service points such as hospitals, schools, pharmacies, markets, water points, or public facilities. If no service file is uploaded, the selected bivariate X variable is used as the access/ranking proxy."),
                    fileInput("p11_service_csv", "Optional service-points CSV", accept = c(".csv")),
                    uiOutput("p11_lat_ui"),
                    uiOutput("p11_lon_ui"),
                    selectInput("p11_service_type", "Service type label:", choices = c("Health facilities", "Schools", "Pharmacies", "Markets", "Water points", "Road/transport access", "Digital services", "Other")),
                    radioButtons("p11_desert_rule", "Service desert definition:", choices = c("Highest distance / lowest access tertile" = "tertile", "Above median distance / below median access" = "median"), selected = "tertile"),
                    actionButton("run_p11", "Run Paper 11 Analysis", class = "btn-success btn-block")
                ),
                box(title = "Table 1. APA Accessibility Summary by Region", width = 8, status = "info", tableOutput("p11_summary_table"))
              ),
              fluidRow(
                box(title = "Figure 1. Service Desert / Accessibility Map", width = 7, status = "success", plotOutput("p11_map", height = "560px")),
                box(title = "Figure 2. Accessibility-Outcome Gradient", width = 5, status = "success", plotOutput("p11_gradient_plot", height = "560px"))
              ),
              fluidRow(
                box(title = "Table 2. Service Desert Classification", width = 6, status = "primary", tableOutput("p11_desert_table")),
                box(title = "Table 3. Outcome Model by Accessibility", width = 6, status = "primary", tableOutput("p11_model_table"))
              )
      ),
      tabItem(tabName = "paper12_gwr",
              fluidRow(
                box(title = "Paper 12 Controls: Local Determinants / GWR-style Analysis", width = 4, status = "primary", solidHeader = TRUE,
                    helpText("This module estimates local neighborhood-specific slopes. It is designed as a robust GWR-style approximation that works even when full GWR packages are not installed."),
                    selectInput("p12_predictor", "Predictor for local coefficient mapping:", choices = NULL),
                    sliderInput("p12_k", "Nearest neighbors for local models:", min = 3, max = 10, value = 5, step = 1),
                    actionButton("run_p12", "Run Paper 12 Analysis", class = "btn-success btn-block")
                ),
                box(title = "Table 1. Global vs Local Model Summary", width = 8, status = "info", tableOutput("p12_global_table"))
              ),
              fluidRow(
                box(title = "Figure 1. Local Coefficient Map", width = 6, status = "success", plotOutput("p12_coef_map", height = "540px")),
                box(title = "Figure 2. Local R-squared Map", width = 6, status = "success", plotOutput("p12_r2_map", height = "540px"))
              ),
              fluidRow(
                box(title = "Table 2. APA Local Coefficients by Region", width = 12, status = "primary", tableOutput("p12_local_table"))
              )
      ),
      tabItem(tabName = "paper13_priority",
              fluidRow(
                box(title = "Paper 13 Controls: Spatial Priority Index", width = 4, status = "primary", solidHeader = TRUE,
                    helpText("Select indicators where higher values indicate higher need or vulnerability. The outcome is included by default."),
                    selectInput("p13_indicators", "Additional indicators for priority score:", choices = NULL, multiple = TRUE),
                    checkboxInput("p13_include_outcome", "Include mapped outcome in priority index", TRUE),
                    actionButton("run_p13", "Run Paper 13 Analysis", class = "btn-success btn-block")
                ),
                box(title = "Table 1. APA Spatial Priority Ranking", width = 8, status = "info", tableOutput("p13_priority_table"))
              ),
              fluidRow(
                box(title = "Figure 1. Priority Index Map", width = 7, status = "success", plotOutput("p13_priority_map", height = "560px")),
                box(title = "Figure 2. Indicator Contribution Summary", width = 5, status = "success", plotOutput("p13_contribution_plot", height = "560px"))
              ),
              fluidRow(
                box(title = "Table 2. Indicator Weights and Contributions", width = 12, status = "primary", tableOutput("p13_indicator_table"))
              )
      ),
      tabItem(tabName = "paper14_scan",
              fluidRow(
                box(title = "Paper 14 Controls: Spatial Scan / Cluster Detection", width = 4, status = "primary", solidHeader = TRUE,
                    helpText("This module combines Global Moran's I, LISA, Getis-Ord Gi*, join-count logic, and contiguous high-burden candidate clusters."),
                    sliderInput("p14_threshold", "High-burden threshold percentile:", min = 50, max = 95, value = 75, step = 5),
                    actionButton("run_p14", "Run Paper 14 Analysis", class = "btn-success btn-block")
                ),
                box(title = "Table 1. APA Global Cluster Diagnostics", width = 8, status = "info", tableOutput("p14_global_table"))
              ),
              fluidRow(
                box(title = "Figure 1. Candidate Spatial Cluster Map", width = 7, status = "success", plotOutput("p14_cluster_map", height = "560px")),
                box(title = "Figure 2. Cluster Burden Profile", width = 5, status = "success", plotOutput("p14_cluster_plot", height = "560px"))
              ),
              fluidRow(
                box(title = "Table 2. APA Local Cluster Classification", width = 12, status = "primary", tableOutput("p14_local_table"))
              )
      ),
      tabItem(tabName = "paper15_equity",
              fluidRow(
                box(title = "Paper 15 Controls: Spatial Equity and Environmental Justice", width = 4, status = "primary", solidHeader = TRUE,
                    helpText("Select an equity stratifier such as sex, wealth, residence, education, displacement, or other vulnerability group."),
                    selectInput("p15_equity_group", "Equity / vulnerability grouping variable:", choices = NULL),
                    uiOutput("p15_vulnerable_ui"),
                    actionButton("run_p15", "Run Paper 15 Analysis", class = "btn-success btn-block")
                ),
                box(title = "Table 1. APA Equity Gap Summary", width = 8, status = "info", tableOutput("p15_gap_table"))
              ),
              fluidRow(
                box(title = "Figure 1. Vulnerable Group Share Map", width = 6, status = "success", plotOutput("p15_share_map", height = "540px")),
                box(title = "Figure 2. Outcome by Equity Group", width = 6, status = "success", plotOutput("p15_group_plot", height = "540px"))
              ),
              fluidRow(
                box(title = "Table 2. APA Regional Equity Burden Table", width = 12, status = "primary", tableOutput("p15_region_table"))
              )
      )

      # --- END NEW TAB ---
    )
  )
)

server <- function(input, output, session) {
  raw_data <- reactive({
    req(input$data_file)
    ext <- tools::file_ext(input$data_file$name)
    if(ext == "dta") read_dta(input$data_file$datapath) else
      if(ext == "sav") read_sav(input$data_file$datapath) else read_csv(input$data_file$datapath)
  })

  observe({
    df <- raw_data(); vars <- names(df)
    updateSelectInput(session, "v_dv", choices = vars)
    updateSelectInput(session, "v_reg", choices = vars)
    updateSelectInput(session, "v_cluster", choices = vars)
    updateSelectInput(session, "v_wgt", choices = vars)
    updateSelectInput(session, "v_strata", choices = vars)
    updateSelectInput(session, "v_indiv_cont", choices = vars)
    updateSelectInput(session, "v_indiv_fact", choices = vars)
    updateSelectInput(session, "v_comm_cont", choices = vars)
    updateSelectInput(session, "v_comm_fact", choices = vars)
    updateSelectInput(session, "v_wealth", choices = vars)
    updateSelectInput(session, "v_time", choices = c("", vars)) # Added "" for optionality
    updateSelectInput(session, "v_risk_group", choices = c("", vars))
    updateSelectInput(session, "joint_dv1", choices = vars)
    updateSelectInput(session, "joint_dv2", choices = vars)
    updateSelectInput(session, "jd_dv1", choices = vars)
    updateSelectInput(session, "jd_dv2", choices = vars)
    updateSelectInput(session, "jd_group", choices = vars)
    updateSelectInput(session, "p12_predictor", choices = vars)
    updateSelectInput(session, "p13_indicators", choices = vars)
    updateSelectInput(session, "p15_equity_group", choices = vars)
  })

  processed_df <- reactive({
    req(input$v_dv)
    df <- raw_data()

    # FORCE conversion of Region to Labels (e.g., 11 -> "Awdal")
    if (!is.null(input$v_reg)) {
      # as_factor() gets the Stata labels, as.character() ensures it's a string
      df[[input$v_reg]] <- as.character(haven::as_factor(df[[input$v_reg]]))
    }

    # Now safely remove labels from other variables for calculation
    df <- df %>% mutate(across(everything(), ~if(inherits(., "haven_labelled")) zap_labels(.) else .))

    y <- as.numeric(as.character(df[[input$v_dv]]))

    if(input$outcome_type == "Binary") {
      if (all(y %in% c(1, 2), na.rm = TRUE)) { y <- y - 1 }
    } else if(input$outcome_type == "Fractional") {
      y <- (y * (length(y) - 1) + 0.5) / length(y)
    } else if(input$outcome_type %in% c("Ordinal", "Multinomial")) {
      y <- as.numeric(as.factor(y))
    }

    df$dv_final <- y
    for(v in c(input$v_indiv_fact, input$v_comm_fact)) { df[[v]] <- as_factor(df[[v]]) }
    return(df)
  })

  # 1. Create formal Survey Design Object
  design_obj <- reactive({
    req(processed_df(), input$v_cluster, input$v_strata, input$v_wgt)
    df <- processed_df() %>% drop_na(dv_final, !!sym(input$v_cluster), !!sym(input$v_strata))
    df$wgt_final <- if(input$is_dhs_wgt) df[[input$v_wgt]]/1000000 else as.numeric(df[[input$v_wgt]])

    survey::svydesign(
      ids = as.formula(paste0("~", input$v_cluster)),
      strata = as.formula(paste0("~", input$v_strata)),
      weights = ~wgt_final,
      data = df,
      nest = TRUE
    )
  })
  #############################################################################################################
  # 2. SAE Engine (Fay-Herriot)
  # --- UPDATED SAE ENGINE (Fixes the 'round' error) ---
  sae_res <- eventReactive(input$run_sae, {
    req(adal_spatial(), input$v_wealth)

    df <- adal_spatial() %>%
      st_drop_geometry() %>%
      filter(!is.na(outcome), !is.na(bi_var), !is.na(sampling_var))

    # Fit Fay-Herriot Model
    fh_mod <- sae::mseFH(outcome ~ bi_var, vardir = sampling_var, data = df)

    # Helper to safely handle NULL/Non-numeric values
    safe_round <- function(x, d = 2) {
      if (is.null(x) || !is.numeric(x) || is.na(x)) return("N/A")
      return(as.character(round(x, d)))
    }

    # Capture Diagnostics for the Paper
    diag_df <- data.frame(
      Metric = c("Log-Likelihood", "AIC", "BIC", "Convergence"),
      Value = c(safe_round(fh_mod$model$fit$logL),
                safe_round(fh_mod$model$fit$AIC),
                safe_round(fh_mod$model$fit$BIC),
                "Successful")
    )

    df$eblup <- fh_mod$est$eblup
    df$mse_eblup <- fh_mod$mse
    df$cv_direct <- (sqrt(df$sampling_var) / df$outcome) * 100
    df$cv_eblup <- (sqrt(df$mse_eblup) / df$eblup) * 100

    return(list(data = df, diags = diag_df))
  })

  # Update existing SAE outputs to handle the new list structure
  output$sae_map <- renderPlot({
    req(sae_res())
    df_sae <- sae_res()$data
    shp <- adal_spatial() %>% dplyr::select(Match_Key, geometry)
    map_df <- shp %>% left_join(df_sae, by = "Match_Key")
    ggplot(map_df) + geom_sf(aes(fill = eblup)) + scale_fill_viridis_c(option = "mako") + theme_void()
  })

  output$sae_table <- renderTable({
    req(sae_res())
    sae_res()$data %>% dplyr::select(Region = Match_Key, outcome, eblup, cv_direct, cv_eblup)
  }, digits = 4)

  output$sae_diag_table <- renderTable({ req(sae_res()); sae_res()$diags })

  # --- NEW OUTPUTS FOR THE SAE PAPER ---
  output$sae_diag_table <- renderTable({ req(sae_res()); sae_res()$diags })

  output$sae_shrinkage_plot <- renderPlot({
    req(sae_res())
    df <- sae_res()$data
    ggplot(df) +
      geom_segment(aes(x = outcome, xend = eblup, y = reorder(Match_Key, outcome), yend = Match_Key),
                   arrow = arrow(length = unit(0.2, "cm")), color = "gray") +
      geom_point(aes(x = outcome, y = Match_Key, color = "Direct Survey"), size = 3) +
      geom_point(aes(x = eblup, y = Match_Key, color = "SAE EBLUP"), size = 3) +
      scale_color_manual(values = c("Direct Survey" = "red", "SAE EBLUP" = "blue")) +
      labs(title = "Shrinkage Effect: How the Model Corrects Survey Estimates",
           subtitle = "Arrows show the shift from noisy direct data to stabilized model estimates",
           x = "Prevalence Estimate", y = "") +
      theme_minimal()
  })
  ###################################################################################################

  # --- FIXED: Spatial Join Logic with Smart Matching (Used for ESDA/Bivariate/ME) ---
  adal_spatial <- eventReactive(input$process_data, {
    req(input$shp_file, design_obj())
    temp_dir <- file.path(tempdir(), as.numeric(Sys.time()))
    dir.create(temp_dir)
    purrr::walk2(input$shp_file$datapath, input$shp_file$name, ~file.copy(.x, file.path(temp_dir, .y)))
    shp <- st_read(list.files(temp_dir, pattern = "\\.shp$", full.names = TRUE)[1]) %>% st_make_valid() %>% st_transform(32638)

    # USE SURVEY PACKAGE FOR DESIGN-BASED MEANS AND SE
    design <- design_obj()
    stats_survey <- survey::svyby(
      formula = ~dv_final,
      by = as.formula(paste0("~", input$v_reg)),
      design = design,
      svymean,
      na.rm = TRUE
    )

    # Calculate Bivariate X Mean for SAE covariate
    df_p <- processed_df()
    df_p$wgt_f <- if(input$is_dhs_wgt) df_p[[input$v_wgt]]/1000000 else as.numeric(df_p[[input$v_wgt]])
    stats_bi <- df_p %>%
      group_by(!!sym(input$v_reg)) %>%
      summarise(bi_var = weighted.mean(as.numeric(as.character(!!sym(input$v_wealth))), wgt_f, na.rm=TRUE),
                n_obs = n(),
                cases = sum(dv_final, na.rm=TRUE))

    stats <- stats_survey %>%
      left_join(stats_bi, by = input$v_reg) %>%
      rename(outcome = dv_final, se_design = se) %>%
      mutate(sampling_var = se_design^2,
             Match_Key = str_to_upper(str_trim(as.character(!!sym(input$v_reg)))))

    # Smart Match logic (unchanged)
    shp_char_cols <- names(shp)[sapply(shp, is.character)]
    best_col <- shp_char_cols[1]; max_matches <- 0
    for(col in shp_char_cols) {
      current_matches <- sum(str_to_upper(str_trim(shp[[col]])) %in% stats$Match_Key)
      if(current_matches > max_matches) { max_matches <- current_matches; best_col <- col }
    }
    shp <- shp %>% mutate(Match_Key = str_to_upper(str_trim(as.character(!!sym(best_col)))))
    joined <- shp %>% left_join(stats, by = "Match_Key") %>% filter(!is.na(outcome))
    return(joined)
  })

  # 1. Spatial Econometrics Engine
  # --- Dedicated Spatial Aggregator for Econometrics (All Predictors) ---
  adal_spatial_econ <- reactive({
    req(adal_spatial())
    df_p <- processed_df()

    preds <- unique(c(input$v_indiv_cont, input$v_indiv_fact,
                      input$v_comm_cont, input$v_comm_fact))
    preds <- setdiff(preds, input$v_reg)

    df_p$wgt_f <- if(input$is_dhs_wgt) df_p[[input$v_wgt]]/1000000 else as.numeric(df_p[[input$v_wgt]])

    stats_all <- df_p %>%
      group_by(!!sym(input$v_reg)) %>%
      summarise(across(all_of(preds),
                       ~ weighted.mean(as.numeric(as.character(.x)), wgt_f, na.rm = TRUE)),
                .groups = "drop") %>%
      mutate(Match_Key = str_to_upper(str_trim(as.character(!!sym(input$v_reg)))))

    # FIXED: Added dplyr::select here
    joined <- adal_spatial() %>%
      left_join(stats_all %>% dplyr::select(-!!sym(input$v_reg)), by = "Match_Key")

    return(list(data = joined, pred_names = preds))
  })
  # --- Multivariate Econometrics Engine ---
  econ_res <- eventReactive(input$run_econ, {
    req(adal_spatial_econ())
    res_obj <- adal_spatial_econ()
    df <- res_obj$data %>% filter(!is.na(outcome))
    preds <- res_obj$pred_names

    # Build Formula: outcome ~ pred1 + pred2 + ...
    form <- as.formula(paste("outcome ~", paste(preds, collapse = " + ")))

    # Weights Matrix
    nb <- if(input$econ_nb_type == "Queen's Contiguity") {
      poly2nb(df, queen = TRUE)
    } else if(input$econ_nb_type == "K-Nearest Neighbors (K=3)") {
      knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 3))
    } else {
      knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 5))
    }

    # Add a check to ensure no islands exist
    if(any(card(nb) == 0)) {
      nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))
    }
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
    # 1. Lagrange Multiplier Tests
    ols_mod <- lm(form, data = df)
    lm_tests <- lm.RStests(ols_mod, lw, test=c("LMerr", "LMlag", "RLMerr", "RLMlag"))
    # 2. Fit Spatial Model
    if(input$econ_model_type == "Spatial Lag Model (SAR)") {
      mod <- lagsarlm(form, data = df, listw = lw, zero.policy = TRUE)
      imp <- impacts(mod, listw = lw)
    } else {
      mod <- errorsarlm(form, data = df, listw = lw, zero.policy = TRUE)
      imp <- NULL
    }

    return(list(model = mod, impacts = imp, lm_tests = lm_tests, pred_names = preds))
  })

  output$econ_lm_tests <- renderPrint({ req(econ_res()); print(econ_res()$lm_tests) })
  output$econ_coef_table <- renderTable({ req(econ_res()); broom::tidy(econ_res()$model) }, striped = TRUE, bordered = TRUE)
  output$econ_impacts_table <- renderTable({
    req(econ_res())
    imp <- econ_res()$impacts
    if(is.null(imp)) return(data.frame(Message = "Impacts only for SAR models."))
    data.frame(Variable = econ_res()$pred_names, Direct = as.numeric(imp$direct), Indirect = as.numeric(imp$indirect), Total = as.numeric(imp$total))
  }, striped = TRUE, bordered = TRUE)

  econ_weights_obj <- reactive({
    req(adal_spatial_econ())
    df <- adal_spatial_econ()$data %>% dplyr::filter(!is.na(outcome))
    safe_spatial_weights(df, input$econ_nb_type)
  })

  econ_fit_table_data <- reactive({
    req(econ_res(), econ_weights_obj())
    mod <- econ_res()$model
    lw <- econ_weights_obj()$lw
    resids <- tryCatch(residuals(mod), error = function(e) NULL)
    moran_res <- tryCatch(spdep::moran.test(resids, lw, zero.policy = TRUE), error = function(e) NULL)
    coef_tbl <- tryCatch(broom::tidy(mod), error = function(e) data.frame())
    spatial_par <- coef_tbl %>% dplyr::filter(term %in% c("rho", "lambda"))
    data.frame(
      Metric = c("Model type", "Spatial weights", "Number of regions", "AIC", "BIC", "Log-likelihood", "Spatial parameter", "Spatial parameter p-value", "Residual Moran's I", "Residual Moran p-value"),
      Value = c(
        input$econ_model_type,
        input$econ_nb_type,
        nrow(adal_spatial_econ()$data),
        round(AIC(mod), 3),
        round(BIC(mod), 3),
        round(as.numeric(logLik(mod)), 3),
        ifelse(nrow(spatial_par) > 0, paste0(spatial_par$term[1], " = ", round(spatial_par$estimate[1], 4)), "Not reported"),
        ifelse(nrow(spatial_par) > 0 && "p.value" %in% names(spatial_par), format_p_apa(spatial_par$p.value[1]), "-"),
        ifelse(is.null(moran_res), "-", round(as.numeric(moran_res$estimate[1]), 4)),
        ifelse(is.null(moran_res), "-", format_p_apa(moran_res$p.value))
      )
    )
  })

  econ_weights_table_data <- reactive({
    req(econ_weights_obj())
    nb <- econ_weights_obj()$nb
    cards <- spdep::card(nb)
    data.frame(
      Metric = c("Regions", "Total neighbor links", "Minimum neighbors", "Mean neighbors", "Maximum neighbors", "Isolated regions"),
      Value = c(length(nb), sum(cards), min(cards), round(mean(cards), 2), max(cards), sum(cards == 0))
    )
  })

  econ_resid_table_data <- reactive({
    req(econ_res(), adal_spatial_econ())
    df <- adal_spatial_econ()$data %>% sf::st_drop_geometry()
    mod <- econ_res()$model
    fitted_vals <- tryCatch(as.numeric(fitted(mod)), error = function(e) rep(NA_real_, nrow(df)))
    resid_vals <- tryCatch(as.numeric(residuals(mod)), error = function(e) rep(NA_real_, nrow(df)))
    df %>%
      dplyr::mutate(Fitted = fitted_vals, Residual = resid_vals, Abs_Residual = abs(Residual)) %>%
      dplyr::arrange(desc(Abs_Residual)) %>%
      dplyr::select(Region = Match_Key, Observed = outcome, Fitted, Residual, Abs_Residual, n_obs, cases) %>%
      dplyr::mutate(across(where(is.numeric), ~round(.x, 4)))
  })

  output$econ_fit_table <- renderTable({ econ_fit_table_data() }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m")
  output$econ_weights_table <- renderTable({ econ_weights_table_data() }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m")
  output$econ_resid_table <- renderTable({ econ_resid_table_data() }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m")


  # --- NEW: Spatiotemporal Data Aggregation ---
  adal_spatial and temporal <- reactive({
    req(input$v_reg, processed_df())
    df_proc <- processed_df()

    if (is.null(input$v_time) || input$v_time == "") return(NULL)

    stats_time <- df_proc %>%
      mutate(
        dv_num = dv_final,
        wgt = if(input$is_dhs_wgt) !!sym(input$v_wgt)/1000000 else as.numeric(!!sym(input$v_wgt))
      ) %>%
      group_by(!!sym(input$v_reg), time_var = as.factor(!!sym(input$v_time))) %>%
      summarise(outcome = weighted.mean(dv_num, wgt, na.rm = TRUE), .groups = 'drop') %>%
      mutate(Match_Key = str_to_upper(str_trim(as.character(!!sym(input$v_reg)))))

    # FIXED: Added dplyr::select here
    shp_base <- adal_spatial() %>%
      dplyr::select(Match_Key, geometry) %>%
      distinct()

    st_sf <- stats_time %>%
      left_join(shp_base, by = "Match_Key") %>%
      st_as_sf()

    return(st_sf)
  })

  output$sparsity_alert <- renderUI({
    df <- adal_spatial()
    low_info <- sum(df$n_obs < 5, na.rm = TRUE)
    percent_low <- (low_info / nrow(df)) * 100
    if(percent_low > 25) {
      callout(title = "Data Sparsity Warning", status = "danger", paste0(round(percent_low, 1), "% of regions have low data."))
    } else {
      callout(title = "Data Density Check", status = "success", "Data distribution is sufficient.")
    }
  })
  # ORIGINAL SPATIAL LOGIC
  generate_esda_plot <- function(df, method) {
    nb <- poly2nb(df, queen = TRUE)
    if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

    if(method == "1. Connectivity Graph") {
      plot(st_geometry(df), border="grey")
      plot(nb, st_coordinates(st_centroid(df)), add=TRUE, col="#006400", lwd=2)
    } else if(method == "2. Basic Choropleth") {
      ggplot(df) + geom_sf(aes(fill=outcome), color="black", size=0.2) +
        scale_fill_distiller(palette = "RdPu", direction = 1, labels = percent) + theme_void()
    } else if(method == "3. EBS Smoothing") {
      eb_res <- EBlocal(ri = df$cases, ni = df$n_obs, nb = nb, zero.policy = TRUE)
      df$eb_rate <- eb_res$est

      # Plot 1: Raw Map with Labels
      p1 <- ggplot(df) +
        geom_sf(aes(fill=outcome)) +
        geom_sf_text(aes(label = Match_Key), size = 2.5, color = "black", check_overlap = TRUE) +
        scale_fill_distiller(palette="RdPu", direction=1) +
        labs(title="Raw Rates") + theme_void()

      # Plot 2: Smoothed Map with Labels
      p2 <- ggplot(df) +
        geom_sf(aes(fill=eb_rate)) +
        geom_sf_text(aes(label = Match_Key), size = 2.5, color = "black", check_overlap = TRUE) +
        scale_fill_distiller(palette="RdPu", direction=1) +
        labs(title="EBS Smoothed") + theme_void()

      plot_grid(p1, p2)
    } else if(method == "4. Quantile Map") {
      breaks_qt <- classIntervals(df$outcome, n = 5, style = "quantile")
      df$quantile_cat <- cut(df$outcome, breaks_qt$brks, include.lowest = TRUE)
      ggplot(df) + geom_sf(aes(fill = quantile_cat), color = "black", size = 0.2) +
        scale_fill_brewer(palette = "Blues", name = "Quintiles") + labs(title = "4. Quantile Map (Ranking)") + theme_void()
    } else if(method == "5. Global Moran's I") {
      gmoran <- moran.test(df$outcome, lw, zero.policy = TRUE)
      z_score <- gmoran$statistic
      p_value <- gmoran$p.value
      moran_i <- gmoran$estimate[1]
      df_norm <- data.frame(x = seq(-4, 4, length.out = 1000)) %>% mutate(y = dnorm(x))
      ggplot(df_norm, aes(x, y)) + geom_line(size = 1) + geom_area(data = subset(df_norm, x >= z_score), fill = "red", alpha = 0.5) +
        geom_vline(xintercept = z_score, color = "red", linetype="dashed") +
        annotate("text", x = 2.5, y = 0.3, label = paste("p-value:", format.pval(p_value)), color="red", fontface="bold") +
        labs(title = "5. Global Moran's I (Systemic Test)", subtitle = paste("Moran's I:", round(moran_i, 3), "| Z-Score:", round(z_score, 2)), x = "Z-Score", y = "Density") + theme_minimal()
    } else if(method == "6. Moran Scatterplot") {
      df$lag_outcome <- lag.listw(lw, df$outcome, zero.policy = TRUE)
      m_prev <- mean(df$outcome)
      df$quadrant <- case_when(df$outcome >= m_prev & df$lag_outcome >= m_prev ~ "High-High", df$outcome <= m_prev & df$lag_outcome <= m_prev ~ "Low-Low", df$outcome >= m_prev & df$lag_outcome <= m_prev ~ "High-Low", df$outcome <= m_prev & df$lag_outcome >= m_prev ~ "Low-High")
      ggplot(df, aes(x = outcome, y = lag_outcome, color = quadrant)) + geom_vline(xintercept = m_prev, linetype = "dashed") + geom_hline(yintercept = mean(df$lag_outcome), linetype = "dashed") + geom_point(size = 3) + geom_text(aes(label = Match_Key), vjust = -0.5, size = 3, show.legend = FALSE) + geom_smooth(method = "lm", se = FALSE, color = "black") + scale_color_manual(values = c("High-High"="red", "Low-Low"="blue", "High-Low"="pink", "Low-High"="lightblue")) + labs(title = "6. Moran Scatterplot (Association)") + theme_minimal()
    } else if(method == "7. LISA Cluster Map") {
      loc_m <- localmoran(df$outcome, lw, zero.policy = TRUE)
      df$Pr <- loc_m[,5]; df$lag_outcome <- lag.listw(lw, df$outcome, zero.policy = TRUE)
      m_prev <- mean(df$outcome)

      # --- ENHANCED LISA CLUSTER MAP (Coloring by Quadrant and Significance) ---
      df$lisa_quadrant <- case_when(
        df$Pr <= 0.05 & df$outcome >= m_prev & df$lag_outcome >= m_prev ~ "High-High",
        df$Pr <= 0.05 & df$outcome <= m_prev & df$lag_outcome <= m_prev ~ "Low-Low",
        df$Pr <= 0.05 & df$outcome >= m_prev & df$lag_outcome <= m_prev ~ "High-Low",
        df$Pr <= 0.05 & df$outcome <= m_prev & df$lag_outcome >= m_prev ~ "Low-High",
        TRUE ~ "Not Significant"
      )

      lisa_colors <- c("High-High"="red", "Low-Low"="blue", "High-Low"="pink", "Low-High"="lightblue", "Not Significant"="lightgray")

      ggplot(df) + geom_sf(aes(fill=lisa_quadrant), color="black") +
        scale_fill_manual(values=lisa_colors, name="LISA Cluster Type") +
        labs(title = "7. LISA Cluster Map (Local Spatial Association)") +
        theme_void()

    } else if(method == "8. Getis-Ord Gi*") {
      gi_stat <- localG(df$outcome, lw, zero.policy = TRUE)
      df$gi_z <- as.numeric(gi_stat)
      ggplot(df) + geom_sf(aes(fill=gi_z)) + scale_fill_gradient2(low="blue", mid="white", high="red") + theme_void()
    } else if(method == "9. Join Count Analysis") {
      global_mean <- mean(df$outcome, na.rm = TRUE)
      df$Coverage_Class <- factor(ifelse(df$outcome > global_mean, "High", "Low"))
      ggplot(data = df) + geom_sf(aes(fill = Coverage_Class), color = "black", size = 0.2) + scale_fill_manual(values = c("High" = "forestgreen", "Low" = "darkorange")) + labs(title = "9. Join Count Binary Map (Categorical Clustering)") + theme_void()
    } else if(method == "10. Hinge Map") {
      cuts <- quantile(df$outcome, c(0, 0.25, 0.5, 0.75, 1)); iqr <- diff(cuts[c(2, 4)])
      df$Hinge <- cut(df$outcome, breaks = c(-Inf, cuts[2]-1.5*iqr, cuts[2], cuts[3], cuts[4], cuts[4]+1.5*iqr, Inf), labels = c("L-Outlier", "<25%", "25-50%", "50-75%", ">75%", "U-Outlier"))
      ggplot(df) + geom_sf(aes(fill=Hinge)) + scale_fill_brewer(palette="RdBu") + theme_void()
    } else if(method == "11. Spatial Correlogram") {
      # Dynamic Lag Calculation:
      n_regions <- nrow(df)
      max_lag <- ifelse(n_regions < 8, floor(n_regions / 2), 4)

      # Safety check: If data is too small (e.g., < 3 regions), show error
      if (max_lag < 1) {
        ggplot() + annotate("text", x = 1, y = 1, label = "Not enough regions for Correlogram", size = 6) + theme_void()
      } else {
        tryCatch({
          cor_res <- sp.correlogram(nb, df$outcome, order = max_lag, method = "I", style = "W", zero.policy = TRUE)

          actual_lags <- nrow(cor_res$res)

          cor_df <- data.frame(
            Lag = 1:actual_lags,
            Morans_I = cor_res$res[, 1],
            P_Value = 2 * (1 - pnorm(abs(cor_res$res[, 1] - cor_res$res[, 2]) / cor_res$res[, 3]))
          ) %>% mutate(Significance = ifelse(P_Value < 0.05, "Significant", "Not Significant"))

          ggplot(cor_df, aes(x = Lag, y = Morans_I)) +
            geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
            geom_line() +
            geom_point(aes(color = Significance), size = 4) +
            scale_color_manual(values = c("Significant" = "red", "Not Significant" = "blue")) +
            scale_x_continuous(breaks = 1:actual_lags) + # Ensure integer ticks
            labs(title = "11. Spatial Correlogram", y = "Moran's I") +
            theme_minimal()
        }, error = function(e) {
          ggplot() + annotate("text", x = 1, y = 1, label = "Data structure cannot support Correlogram", size = 5) + theme_void()
        })
      }} else if(method == "12. Regime Density") {
        gi_stat <- localG(df$outcome, lw, zero.policy = TRUE)
        df$Regime <- ifelse(as.numeric(gi_stat) > 1.65, "Hot Spot", ifelse(as.numeric(gi_stat) < -1.65, "Cold Spot", "Neutral"))
        ggplot(df, aes(x=outcome, fill=Regime)) + geom_density(alpha=0.5) + theme_minimal()
      } else if(method == "13. Spatial Lag Map") {
        df$lag_outcome <- lag.listw(lw, df$outcome, zero.policy = TRUE)
        ggplot(df) + geom_sf(aes(fill = lag_outcome), color = "black", size = 0.2) + scale_fill_distiller(palette = "YlGnBu", direction = 1, labels = percent) + labs(title = "13. Spatial Lag Map (Neighborhood Context)", subtitle = "Average Rate of Neighbors") + theme_void()
      } else if(method == "14. SDG Target Gap") {
        df$gap <- df$outcome - 0
        ggplot(df) + geom_sf(aes(fill=gap)) + scale_fill_viridis_c(option="magma", direction=-1) + theme_void()
      } else if(method == "15. Local Geary's C") {
        geary_loc <- localC(df$outcome, lw); df$geary_c <- geary_loc
        ggplot(df) + geom_sf(aes(fill = geary_c), color = "black", size = 0.2) + scale_fill_viridis_c(option = "magma", direction = -1) + labs(title = "15. Local Geary's C (Edge Detection)", subtitle = "Darker = High Dissimilarity") + theme_void()
      } else if(method == "16. Conditional Plot") {
        df$lag_bin <- cut(lag.listw(lw, df$outcome), breaks=3, labels=c("Low Nb", "Med Nb", "High Nb"))
        ggplot(df, aes(x=lag_bin, y=outcome, fill=lag_bin)) + geom_boxplot() + geom_jitter(width=0.1) + theme_minimal()
      } else if(method == "17. Global Geary's C") {
        g_geary <- geary.test(df$outcome, lw, zero.policy = TRUE)
        z_score <- as.numeric(g_geary$statistic)
        p_value <- g_geary$p.value
        geary_c <- g_geary$estimate[1]
        df_norm <- data.frame(x = seq(-4, 4, length.out = 1000)) %>% mutate(y = dnorm(x))
        ggplot(df_norm, aes(x, y)) + geom_line(size = 1) +
          geom_area(data = subset(df_norm, x <= z_score), fill = "#006400", alpha = 0.5) +
          geom_vline(xintercept = z_score, color = "red", linetype="dashed", size = 1) +
          annotate("text", x = 2.5, y = 0.3, label = paste("p-value:", format.pval(p_value)), color="red", fontface="bold") +
          labs(title = "17. Global Geary's C (Systemic Test)",
               subtitle = paste("Geary's C:", round(geary_c, 3), "| Z-Score:", round(z_score, 2)),
               x = "Z-Score", y = "Density") + theme_minimal()

      } else {
        ggplot(df) + geom_sf(aes(fill=outcome)) + scale_fill_viridis_c() + theme_void()
      }
  }

  output$diag_plot <- renderPlot({
    df <- adal_spatial()
    generate_esda_plot(df, input$esda_method)
  })

  output$bi_moran_test <- renderPrint({
    df <- adal_spatial()
    req(input$v_wealth)

    # 1. Define Spatial Weights
    nb <- poly2nb(df, queen = TRUE)
    if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

    # Check for NA values in the variables used for the test
    if (any(is.na(df$outcome)) || any(is.na(df$bi_var))) {
      cat("Bivariate Moran's I Test Skipped: Outcome or Bivariate variable contains missing values (NA).")
      return(NULL)
    }

    # --- CORRECTED CONDITION CHECK ---
    # Check if we have enough observations AND if the spatial weights list is valid
    is_valid_spatial_structure <- nrow(df) > 1 && length(lw$neighbours) > 1 && sum(card(nb)) > 0

    if (is_valid_spatial_structure) {

      # 2. Bivariate Moran's I (Y vs Lagged X) using Moran Monte Carlo (the robust method)

      # Standardize variables for Bivariate I calculation
      z_y <- scale(df$outcome)[, 1]
      z_x <- scale(df$bi_var)[, 1]

      # The Bivariate Moran's I statistic is calculated on the product of the standardized variables
      # and tested via permutation (Moran Monte Carlo).
      product_var <- z_y * z_x

      # Run Moran Monte Carlo test (99 simulations)
      bi_moran_mc <- moran.mc(product_var, lw, nsim = 99, zero.policy = TRUE)

      # Extract results
      I_B_observed <- bi_moran_mc$statistic
      p_value <- bi_moran_mc$p.value

      cat("Bivariate Moran's I (Outcome vs. Lagged X) - Permutation Test (999 simulations):\n")
      cat(paste("Observed Bivariate I:", round(I_B_observed, 4), "\n"))
      cat(paste("P-value (Simulated):", format.pval(p_value), "\n"))
      cat(paste("Alternative hypothesis: greater\n"))

      cat("\nInterpretation: Measures if the outcome in a region is correlated with the X variable (e.g., Wealth) in neighboring regions. A significant positive I indicates spatial covariance between the two variables.")

    } else {
      cat("Bivariate Moran's I Test Skipped: Requires at least two connected regions (e.g., check shapefile connectivity or data density).")
    }
    # --- END CORRECTED CONDITION CHECK ---
  })

  output$bi_plot <- renderPlot({
    req(input$run_bi)
    df <- adal_spatial()
    bi_data <- bi_class(df, x = bi_var, y = outcome, style = "quantile", dim = 3)
    map <- ggplot(bi_data) + geom_sf(aes(fill = bi_class), show.legend = FALSE) +
      bi_scale_fill(pal = "GrPink", dim = 3) + theme_void()
    legend <- bi_legend(pal = "GrPink", dim = 3, xlab = "Higher X ", ylab = "Higher Outcome ")
    ggdraw() + draw_plot(map, 0, 0, 1, 1) + draw_plot(legend, 0.05, 0.05, 0.25, 0.25)
  })


  output$krig_plot <- renderPlot({
    df <- adal_spatial()
    # Add st_jitter to prevent singular matrix errors from overlapping centroids
    pts <- st_centroid(df) %>% st_jitter(amount = 0.01)

    grid <- st_make_grid(df, n = c(50, 50), what = "centers") %>% st_as_sf() %>% st_filter(df)
    v_fit <- variogram(outcome ~ 1, pts)
    v_mod <- fit.variogram(v_fit, model = vgm(psill=var(pts$outcome, na.rm=T), "Exp", range=500000, nugget=0.01))
    krig_res <- krige(outcome ~ 1, pts, grid, model = v_mod)
    ggplot(krig_res) + geom_sf(aes(color=var1.pred), size=1) + scale_color_viridis_c() + theme_void()
  })
  output$resid_plot <- renderPlot({
    df <- adal_spatial(); nb <- poly2nb(df, queen = TRUE); lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
    m_sar <- lagsarlm(outcome ~ 1, data = df, lw, zero.policy = TRUE)
    df$resids <- residuals(m_sar)
    ggplot(df) + geom_sf(aes(fill=resids)) + scale_fill_gradient2() + theme_void()
  })


  # --- ADDED FOR FINAL REPUBLISH VERSION: Bivariate screening table in Advanced Spatial tab ---
  bivariate_screening_table <- reactive({
    req(processed_df())
    calc_bivariate_screening_table(
      df = processed_df(),
      dv_col = "dv_final",
      cont_vars = c(input$v_indiv_cont, input$v_comm_cont),
      cat_vars = c(input$v_indiv_fact, input$v_comm_fact)
    )
  })

  output$advanced_biv_table <- renderTable({
    bivariate_screening_table()
  }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m", width = "100%")
  ###############################################################################################################
  # 1. Calculate Inequality Indices
  ineq_stats <- reactive({
    req(adal_spatial())
    df <- adal_spatial() %>% filter(!is.na(outcome))

    g_val <- ineq(df$outcome, type = "Gini")
    t_val <- ineq(df$outcome, type = "Theil")

    return(list(gini = g_val, theil = t_val))
  })

  # 2. Render InfoBoxes
  output$gini_box <- renderInfoBox({
    infoBox("Gini Index", round(ineq_stats()$gini, 3), icon = icon("divide"), color = "danger", fill = TRUE) # Changed 'red' to 'danger'
  })

  output$theil_box <- renderInfoBox({
    infoBox("Theil Index", round(ineq_stats()$theil, 3), icon = icon("chart-line"), color = "orange", fill = TRUE)
  })

  # 3. Lorenz Curve Plot
  output$lorenz_plot <- renderPlot({
    req(adal_spatial())
    df <- adal_spatial() %>% filter(!is.na(outcome))

    # Plotting the Lorenz Curve
    plot(Lc(df$outcome), col = "red", lwd = 2, main = "Geographic Lorenz Curve")
  })
  # Add to ineq_stats reactive
  conc_index <- reactive({
    req(adal_spatial(), input$v_wealth)
    df <- adal_spatial() %>% st_drop_geometry() %>% filter(!is.na(outcome), !is.na(bi_var))
    # bi_var is your wealth variable aggregated to region
    # We sort by wealth to calculate the Concentration Index
    df <- df[order(df$bi_var), ]
    # Simple Concentration Index calculation
    mu <- mean(df$outcome)
    n <- nrow(df)
    index <- (2 / (n * mu)) * sum(df$outcome * (1:n)) - 1 - (1/n)
    return(index)
  })

  output$conc_box <- renderInfoBox({
    infoBox("Concentration Index", round(conc_index(), 4),
            icon = icon("hand-holding-usd"), color = "purple", fill = TRUE)
  })


  # --- ADDED FOR FINAL REPUBLISH VERSION: Detailed APA inequality tables ---
  inequality_apa_tables <- reactive({
    calc_inequality_apa_tables(adal_spatial(), wealth_available = !is.null(input$v_wealth) && input$v_wealth != "")
  })

  output$ineq_apa_table <- renderTable({
    req(inequality_apa_tables())
    inequality_apa_tables()$summary %>%
      dplyr::select(Metric, Estimate = APA_Formatted, Interpretation)
  }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m", width = "100%")

  output$ineq_region_detail_table <- renderTable({
    req(inequality_apa_tables())
    inequality_apa_tables()$regional %>%
      dplyr::mutate(
        Outcome = round(Outcome, 4),
        Outcome_Share_Percent = round(Outcome_Share_Percent, 2),
        Cumulative_Share_Percent = round(Cumulative_Share_Percent, 2),
        Gap_From_Mean = round(Gap_From_Mean, 4),
        Relative_Gap_Percent = round(Relative_Gap_Percent, 2)
      )
  }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m", width = "100%")
  # --- SPATIOTEMPORAL OUTPUTS ---

  output$st_moran_trend_ui <- renderUI({
    df_st <- adal_spatial and temporal()

    if (is.null(df_st) || length(unique(df_st$time_var)) < 2) {
      bs4Dash::callout(
        title = "Spatiotemporal Analysis Disabled",
        status = "warning",
        width = 12,
        "This module requires at least two survey waves (years). Please ensure a Time Variable is selected in the Setup tab."
      )
    } else {
      plotOutput("st_moran_trend", height = "400px")
    }
  })

  st_moran_trend_results <- reactive({
    df_st <- adal_spatial and temporal()
    if (is.null(df_st) || nrow(df_st) == 0) return(NULL)
    time_points_numeric <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique()
    if (length(time_points_numeric) < 2) return(NULL)

    df_st %>%
      group_by(time_var) %>%
      group_modify(~ {
        current_df <- .x
        nb <- poly2nb(current_df, queen = TRUE)
        if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(current_df)), k = 1))
        lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
        if (length(lw$neighbours) > 1) {
          m_test <- moran.test(current_df$outcome, lw, zero.policy = TRUE)
          tibble(
            Moran_I = as.numeric(m_test$estimate[1]),
            Expected_I = as.numeric(m_test$estimate[2]),
            Variance = as.numeric(m_test$estimate[3]),
            Z_value = as.numeric(m_test$statistic),
            P_value_raw = m_test$p.value
          )
        } else {
          tibble(Moran_I = NA_real_, Expected_I = NA_real_, Variance = NA_real_, Z_value = NA_real_, P_value_raw = NA_real_)
        }
      }) %>%
      ungroup() %>%
      mutate(
        P_value = sapply(P_value_raw, format_p_apa),
        Significance = sapply(P_value_raw, sig_label),
        Interpretation = sapply(P_value_raw, interpret_p)
      )
  })

  output$st_moran_trend_table <- renderTable({
    res <- st_moran_trend_results()
    if (is.null(res)) return(data.frame(Message = "Requires at least two time periods."))
    res %>% dplyr::mutate(
      Moran_I = round(Moran_I, 4),
      Expected_I = round(Expected_I, 4),
      Variance = round(Variance, 4),
      Z_value = round(Z_value, 3)
    ) %>% dplyr::select(Time = time_var, Moran_I, Expected_I, Variance, Z_value, P_value, Significance, Interpretation)
  }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m", width = "100%")

  output$st_moran_trend <- renderPlot({
    df_st <- adal_spatial and temporal()
    req(nrow(df_st) > 0)

    # Robustly check unique time points (already checked in UI, but safe here)
    time_points_numeric <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique()
    req(length(time_points_numeric) >= 2)

    # Calculate Moran's I for each time period
    moran_results <- df_st %>%
      group_by(time_var) %>%
      group_modify(~ {
        current_df <- .x
        # Ensure connectivity for the current subset
        nb <- poly2nb(current_df, queen = TRUE)
        # Fallback to KNN if any region is isolated
        if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(current_df)), k = 1))
        lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

        # Calculate Moran's I
        if (length(lw$neighbours) > 1) {
          m_test <- moran.test(current_df$outcome, lw, zero.policy = TRUE)
          tibble(
            Moran_I = m_test$estimate[1],
            P_Value = m_test$p.value
          )
        } else {
          tibble(Moran_I = NA, P_Value = NA)
        }
      }) %>%
      ungroup() %>%
      mutate(Significant = P_Value < 0.05)

    ggplot(moran_results, aes(x = time_var, y = Moran_I, group = 1)) +
      geom_line(color = "#006400", size = 1.2) +
      geom_point(aes(color = Significant), size = 4) +
      geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
      scale_color_manual(values = c("TRUE" = "red", "FALSE" = "gray"), name = "Significant (p<0.05)") +
      labs(title = "Global Spatial Autocorrelation Trend (Moran's I)",
           subtitle = "Tracking the strength of clustering over time.",
           x = input$v_time, y = "Moran's I") +
      theme_minimal(base_size = 14)
  })

  output$st_rate_change_ui <- renderUI({
    df_st <- adal_spatial and temporal()
    time_points_numeric <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique()

    if (length(time_points_numeric) < 2) {
      box(title = "4.3 Outcome Rate of Change (Min to Max Time)", status = "warning", width = 12,
          p("Skipped: Requires multiple time points (survey years) for Rate of Change analysis."))
    } else {
      plotOutput("st_rate_change", height = "400px")
    }
  })

  output$st_rate_change <- renderPlot({
    df_st <- adal_spatial and temporal()

    # --- FIX: Convert factor time_var to numeric for min/max ---
    time_points <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique() %>% sort()

    req(length(time_points) >= 2)

    t_min <- min(time_points)
    t_max <- max(time_points)

    # Filter using the numeric time points, matching against the factor/character column
    df_min <- df_st %>%
      filter(as.numeric(as.character(time_var)) == t_min) %>%
      st_drop_geometry() %>%
      dplyr::select(Match_Key, outcome_min = outcome) # Added dplyr::

    df_max <- df_st %>%
      filter(as.numeric(as.character(time_var)) == t_max) %>%
      st_drop_geometry() %>%
      dplyr::select(Match_Key, outcome_max = outcome) # Added dplyr::

    df_change <- df_min %>%
      full_join(df_max, by = "Match_Key") %>%
      mutate(
        change = outcome_max - outcome_min,
        percent_change = (change / outcome_min) * 100
      )

    map_df <- adal_spatial() %>%
      dplyr::select(Match_Key, geometry) %>% # Added dplyr::
      left_join(df_change, by = "Match_Key")

    ggplot(map_df) +
      geom_sf(aes(fill = percent_change), color = "white", size = 0.2) +
      scale_fill_gradient2(low = "red", mid = "white", high = "blue",
                           name = "Percent Change",
                           labels = function(x) paste0(round(x, 0), "%")) +
      theme_void() +
      labs(title = paste("Percent Change in Outcome Rate:", t_min, "to", t_max),
           subtitle = "Blue = Improvement (Positive Change), Red = Decline (Negative Change)")
  })

  output$st_lisa_transition_ui <- renderUI({
    df_st <- adal_spatial and temporal()
    time_points_numeric <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique()

    if (length(time_points_numeric) < 2) {
      box(title = "4.4 LISA Cluster Transition Map", status = "warning", width = 12,
          p("Skipped: Requires multiple time points (survey years) for LISA Transition analysis."))
    } else {
      plotOutput("st_lisa_transition", height = "400px")
    }
  })

  output$st_lisa_transition <- renderPlot({
    df_st <- adal_spatial and temporal()

    # --- FIX: Convert factor time_var to numeric for min/max ---
    time_points <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique() %>% sort()

    req(length(time_points) >= 2)

    t_min <- min(time_points)
    t_max <- max(time_points)

    # 1. Get base geometry and neighborhood list
    base_shp <- df_st %>% filter(as.numeric(as.character(time_var)) == t_min)
    nb <- poly2nb(base_shp, queen = TRUE)
    if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(base_shp)), k = 1))
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

    # 2. Calculate LISA for Min Time
    df_min <- df_st %>% filter(as.numeric(as.character(time_var)) == t_min)
    loc_m_min <- localmoran(df_min$outcome, lw, zero.policy = TRUE)
    m_prev_min <- mean(df_min$outcome)
    df_min$lisa_min <- case_when(
      loc_m_min[,5] <= 0.05 & df_min$outcome >= m_prev_min ~ "HH",
      loc_m_min[,5] <= 0.05 & df_min$outcome <= m_prev_min ~ "LL",
      TRUE ~ "NS"
    )

    # 3. Calculate LISA for Max Time
    df_max <- df_st %>% filter(as.numeric(as.character(time_var)) == t_max)
    loc_m_max <- localmoran(df_max$outcome, lw, zero.policy = TRUE)
    m_prev_max <- mean(df_max$outcome)
    df_max$lisa_max <- case_when(
      loc_m_max[,5] <= 0.05 & df_max$outcome >= m_prev_max ~ "HH",
      loc_m_max[,5] <= 0.05 & df_max$outcome <= m_prev_max ~ "LL",
      TRUE ~ "NS"
    )

    # 4. Join and create Transition Category
    df_transition <- df_min %>%
      st_drop_geometry() %>%
      dplyr::select(Match_Key, lisa_min) %>% # Added dplyr::
      left_join(df_max %>% st_drop_geometry() %>% dplyr::select(Match_Key, lisa_max), by = "Match_Key") %>%
      mutate(
        Transition = paste0(lisa_min, " -> ", lisa_max),
        Transition = ifelse(lisa_min == "NS" | lisa_max == "NS", "Not Significant", Transition)
      )

    map_df <- adal_spatial() %>%
      dplyr::select(Match_Key, geometry) %>% # Added dplyr::
      left_join(df_transition, by = "Match_Key")
    # Define a color palette for key transitions
    transition_colors <- c(
      "HH -> HH" = "darkred",
      "LL -> LL" = "darkblue",
      "HH -> LL" = "purple",
      "LL -> HH" = "orange",
      "NS -> NS" = "lightgray",
      "Not Significant" = "lightgray"
    )

    ggplot(map_df) +
      geom_sf(aes(fill = Transition), color = "black", size = 0.1) +
      scale_fill_manual(values = transition_colors, na.value = "white") +
      theme_void() +
      labs(title = paste("LISA Cluster Transition:", t_min, "to", t_max),
           subtitle = "HH = High-High, LL = Low-Low")
  })

  st_rate_change_table_data <- reactive({
    df_st <- adal_spatial and temporal()
    if (is.null(df_st) || nrow(df_st) == 0) return(data.frame(Message = "Spatiotemporal data not available."))
    time_points <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique() %>% sort()
    if (length(time_points) < 2) return(data.frame(Message = "Requires at least two time periods."))
    t_min <- min(time_points); t_max <- max(time_points)
    df_min <- df_st %>% filter(as.numeric(as.character(time_var)) == t_min) %>% st_drop_geometry() %>% dplyr::select(Match_Key, outcome_min = outcome)
    df_max <- df_st %>% filter(as.numeric(as.character(time_var)) == t_max) %>% st_drop_geometry() %>% dplyr::select(Match_Key, outcome_max = outcome)
    df_min %>% full_join(df_max, by = "Match_Key") %>%
      mutate(
        Absolute_Change = outcome_max - outcome_min,
        Percent_Change = ifelse(outcome_min == 0, NA_real_, (Absolute_Change / outcome_min) * 100),
        Direction = case_when(Absolute_Change > 0 ~ "Increase", Absolute_Change < 0 ~ "Decrease", TRUE ~ "No change"),
        Rank_By_Absolute_Change = rank(-abs(Absolute_Change), ties.method = "average")
      ) %>%
      arrange(Rank_By_Absolute_Change) %>%
      mutate(across(where(is.numeric), ~round(.x, 4))) %>%
      dplyr::rename(Region = Match_Key, Baseline_Outcome = outcome_min, Endline_Outcome = outcome_max)
  })

  output$st_rate_change_table <- renderTable({ st_rate_change_table_data() }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m")

  st_lisa_transition_table_data <- reactive({
    df_st <- adal_spatial and temporal()
    if (is.null(df_st) || nrow(df_st) == 0) return(data.frame(Message = "Spatiotemporal data not available."))
    time_points <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique() %>% sort()
    if (length(time_points) < 2) return(data.frame(Message = "Requires at least two time periods."))
    t_min <- min(time_points); t_max <- max(time_points)
    base_shp <- df_st %>% filter(as.numeric(as.character(time_var)) == t_min)
    sw <- safe_spatial_weights(base_shp, "Queen's Contiguity")
    lw <- sw$lw
    df_min <- df_st %>% filter(as.numeric(as.character(time_var)) == t_min)
    df_max <- df_st %>% filter(as.numeric(as.character(time_var)) == t_max)
    loc_m_min <- localmoran(df_min$outcome, lw, zero.policy = TRUE)
    loc_m_max <- localmoran(df_max$outcome, lw, zero.policy = TRUE)
    lag_min <- lag.listw(lw, df_min$outcome, zero.policy = TRUE)
    lag_max <- lag.listw(lw, df_max$outcome, zero.policy = TRUE)
    m_min <- mean(df_min$outcome, na.rm = TRUE); m_max <- mean(df_max$outcome, na.rm = TRUE)
    lisa_min <- case_when(
      loc_m_min[,5] <= 0.05 & df_min$outcome >= m_min & lag_min >= m_min ~ "High-High",
      loc_m_min[,5] <= 0.05 & df_min$outcome <= m_min & lag_min <= m_min ~ "Low-Low",
      loc_m_min[,5] <= 0.05 & df_min$outcome >= m_min & lag_min <= m_min ~ "High-Low",
      loc_m_min[,5] <= 0.05 & df_min$outcome <= m_min & lag_min >= m_min ~ "Low-High",
      TRUE ~ "Not significant"
    )
    lisa_max <- case_when(
      loc_m_max[,5] <= 0.05 & df_max$outcome >= m_max & lag_max >= m_max ~ "High-High",
      loc_m_max[,5] <= 0.05 & df_max$outcome <= m_max & lag_max <= m_max ~ "Low-Low",
      loc_m_max[,5] <= 0.05 & df_max$outcome >= m_max & lag_max <= m_max ~ "High-Low",
      loc_m_max[,5] <= 0.05 & df_max$outcome <= m_max & lag_max >= m_max ~ "Low-High",
      TRUE ~ "Not significant"
    )
    data.frame(
      Region = df_min$Match_Key,
      Baseline_LISA = lisa_min,
      Endline_LISA = lisa_max,
      Transition = paste(lisa_min, "to", lisa_max),
      Baseline_P = sapply(loc_m_min[,5], format_p_apa),
      Endline_P = sapply(loc_m_max[,5], format_p_apa),
      stringsAsFactors = FALSE
    )
  })

  output$st_lisa_transition_table <- renderTable({ st_lisa_transition_table_data() }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m")

  st_hotspot_table_data <- reactive({
    res <- tryCatch(st_int_res(), error = function(e) NULL)
    if (is.null(res) || is.null(res$interaction_data)) return(data.frame(Message = "Run the spatial and temporal interaction model first."))
    map_df <- adal_spatial() %>% sf::st_drop_geometry() %>% dplyr::select(Match_Key, outcome, n_obs)
    res$interaction_data %>%
      dplyr::left_join(map_df, by = "Match_Key") %>%
      dplyr::mutate(
        Hotspot_Class = case_when(
          interaction_effect > 0 ~ "Emerging hotspot / worsening faster than average",
          interaction_effect < 0 ~ "Improving faster than average",
          TRUE ~ "No differential trend"
        ),
        Rank = rank(-interaction_effect, ties.method = "average")
      ) %>%
      dplyr::arrange(Rank) %>%
      dplyr::select(Region = Match_Key, Interaction_Effect = interaction_effect, Outcome = outcome, n_obs, Hotspot_Class, Rank) %>%
      dplyr::mutate(across(where(is.numeric), ~round(.x, 4)))
  })

  output$st_hotspot_table <- renderTable({ st_hotspot_table_data() }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m")

  output$st_small_multiples <- renderPlot({
    df_st <- adal_spatial and temporal()
    req(nrow(df_st) > 0)

    # Determine the range for consistent color scaling across all time periods
    max_outcome <- max(df_st$outcome, na.rm = TRUE)
    min_outcome <- min(df_st$outcome, na.rm = TRUE)

    ggplot(df_st) +
      geom_sf(aes(fill = outcome), color = "black", size = 0.1) +
      scale_fill_viridis_c(option = "plasma", labels = percent, limits = c(min_outcome, max_outcome)) +
      facet_wrap(~ time_var, ncol = 3) +
      labs(title = "Outcome Rate Over Time (Small Multiples)",
           fill = "Outcome Rate") +
      theme_void() +
      theme(strip.text = element_text(face = "bold", size = 12),
            legend.position = "bottom")
  })
  ########################################################################################################
  # --- UPDATED: MULTILEVEL 4-STAGE PUBLICATION ENGINE ---

  # --- 1. Multilevel 4-Stage Engine ---
  freq_res <- eventReactive(input$run_freq, {
    req(input$v_dv, input$v_cluster)

    # Define predictor sets
    preds_indiv <- c(input$v_indiv_cont, input$v_indiv_fact)
    preds_comm <- c(input$v_comm_cont, input$v_comm_fact)
    all_preds <- unique(c(preds_indiv, preds_comm))

    # Clean data: remove NAs from DV and Cluster
    df_clean <- processed_df() %>%
      drop_na(dv_final, !!sym(input$v_cluster))

    df_clean$cluster_idx <- as.factor(df_clean[[input$v_cluster]])

    # Map glmmTMB Families based on outcome type
    fam_obj <- switch(input$outcome_type,
                      "Binary"     = binomial(link="logit"),
                      "Count"      = poisson(),
                      "Continuous" = gaussian(),
                      "Fractional" = beta_family(),
                      "Ordinal"    = cumulative(),
                      "Nominal"    = multinomial(),
                      gaussian())

    # Formula builder helper
    build_form_freq <- function(preds) {
      rhs <- "(1 | cluster_idx)"
      if (length(preds) > 0) rhs <- paste(paste(preds, collapse = " + "), "+", rhs)
      as.formula(paste("dv_final ~", rhs))
    }

    withProgress(message = 'Running Hierarchical Models...', value = 0, {
      incProgress(0.2, detail = "Stage 0: Null Model")
      m0 <- glmmTMB(build_form_freq(NULL), data = df_clean, family = fam_obj)

      incProgress(0.2, detail = "Stage 1: Individual Level")
      m1 <- if(length(preds_indiv) > 0) glmmTMB(build_form_freq(preds_indiv), data = df_clean, family = fam_obj) else m0

      incProgress(0.2, detail = "Stage 2: Community Level")
      m2 <- if(length(preds_comm) > 0) glmmTMB(build_form_freq(preds_comm), data = df_clean, family = fam_obj) else m0

      incProgress(0.2, detail = "Stage 3: Full Model")
      m3 <- if(length(all_preds) > 0) glmmTMB(build_form_freq(all_preds), data = df_clean, family = fam_obj) else m0
    })

    list(m0=m0, m1=m1, m2=m2, m3=m3)
  })

  # --- 2. Formal Side-by-Side Publication Table ---
  # --- FIXED: Formal Side-by-Side Publication Table ---
  output$ml_comp_table <- renderTable({
    req(freq_res())
    models <- freq_res()

    # Helper function to tidy and format each model column
    tidy_m_full <- function(m, name) {
      # If the model failed (is an error object), return a dummy dataframe
      if (inherits(m, "try-error") || is.null(m)) {
        return(data.frame(term = character(), placeholder = character()) %>%
                 dplyr::rename(!!name := placeholder))
      }

      # Get fixed effects with 95% Confidence Intervals
      res <- tryCatch({
        broom.mixed::tidy(m, effects = "fixed", conf.int = TRUE)
      }, error = function(e) return(NULL))

      if (is.null(res)) return(data.frame(term = "Error", val = "-") %>% dplyr::rename(!!name := val))

      # For Binary and Count, convert log-odds/log-rates to Ratios (AOR/ARR)
      if(input$outcome_type %in% c("Binary", "Count")) {
        res <- res %>% mutate(across(c(estimate, conf.low, conf.high), exp))
      }

      res %>%
        mutate(
          stars = case_when(
            p.value < 0.001 ~ "***",
            p.value < 0.01  ~ "**",
            p.value < 0.05  ~ "*",
            TRUE            ~ ""
          ),
          # Create the academic string: Estimate (Lower CI, Upper CI) Stars
          formatted = paste0(
            round(estimate, 2),
            " (", round(conf.low, 2), "-", round(conf.high, 2), ")",
            stars
          )
        ) %>%
        # USE DPLYR::RENAME AND DPLYR::SELECT EXPLICITLY TO AVOID CONFLICTS
        dplyr::rename(!!name := formatted) %>%
        dplyr::select(term, !!sym(name))
    }

    # Process all 4 models
    t0 <- tidy_m_full(models$m0, "Model 0 (Null)")
    t1 <- tidy_m_full(models$m1, "Model 1 (Indiv)")
    t2 <- tidy_m_full(models$m2, "Model 2 (Comm)")
    t3 <- tidy_m_full(models$m3, "Model 3 (Full)")

    # Join models side-by-side using purrr::reduce for safety
    final_table <- list(t0, t1, t2, t3) %>%
      purrr::reduce(full_join, by = "term") %>%
      # Remove the Intercept row for a cleaner publication look
      filter(!str_detect(term, "Intercept")) %>%
      dplyr::rename(Predictor = term)

    # Replace NA with "-" for variables not present in specific models
    final_table[is.na(final_table)] <- "-"

    return(final_table)

  }, striped = TRUE, bordered = TRUE, align = 'l', spacing = 'm')

  output$freq_forest_plot <- renderPlot({
    req(freq_res())
    model <- freq_res()$m3
    res_df <- broom.mixed::tidy(model, effects = "fixed", conf.int = TRUE) %>%
      filter(term != "(Intercept)")

    ggplot(res_df, aes(x=estimate, y=term)) +
      geom_point(size=4, color="#006400") +
      geom_errorbarh(aes(xmin=conf.low, xmax=conf.high), color="#FFD700", height=0.2) +
      geom_vline(xintercept = 0, linetype="dashed") +
      theme_minimal() + labs(title="Full Model: Fixed Effects (95% CI)")
  })

  output$freq_metrics_plot <- renderPlot({
    req(freq_res())
    res <- freq_res()

    # 1. Construct the data frame with explicit numeric conversion for LogLik
    metrics_df <- data.frame(
      Model = factor(c("Null", "Indiv", "Comm", "Full"),
                     levels = c("Null", "Indiv", "Comm", "Full")),
      AIC = c(AIC(res$m0), AIC(res$m1), AIC(res$m2), AIC(res$m3)),
      BIC = c(BIC(res$m0), BIC(res$m1), BIC(res$m2), BIC(res$m3)),
      ICC = c(calc_freq_icc(res$m0), calc_freq_icc(res$m1), calc_freq_icc(res$m2), calc_freq_icc(res$m3))
    )

    # Calculate PCV (Proportional Change in Variance)
    pcv_m1 <- calc_pcv(res$m0, res$m1)
    pcv_m2 <- calc_pcv(res$m0, res$m2)
    pcv_m3 <- calc_pcv(res$m0, res$m3)

    pcv_df <- data.frame(
      Model = factor(c("Indiv", "Comm", "Full"), levels = c("Indiv", "Comm", "Full")),
      PCV = c(pcv_m1, pcv_m2, pcv_m3)
    )

    # Helper function for consistent professional styling
    style_metric_plot <- function(df, y_var, title, is_icc = FALSE, is_pcv = FALSE) {
      p <- ggplot(df, aes(x = Model, y = !!sym(y_var), fill = Model)) +
        geom_bar(stat = "identity", width = 0.6, show.legend = FALSE) +
        geom_text(aes(label = round(!!sym(y_var), 3)), vjust = -0.3, size = 3.5, fontface = "bold") +
        scale_fill_viridis_d(option = "mako", begin = 0.3, end = 0.8) +
        labs(title = title, x = NULL, y = NULL) +
        theme_minimal(base_size = 12) +
        theme(
          plot.title = element_text(size = 11, face = "bold", hjust = 0.5),
          panel.grid.major.x = element_blank(),
          axis.text.x = element_text(face = "bold")
        )

      if(is_icc) p <- p + ylim(0, 1.1)
      if(is_pcv) p <- p + labs(y = "PCV (Proportion Explained)") +
          geom_text(aes(label = scales::percent(!!sym(y_var), accuracy = 0.1)), vjust = -0.3, size = 3.5, fontface = "bold")

      return(p)
    }

    # 2. Generate individual plots
    p_aic  <- style_metric_plot(metrics_df, "AIC", "AIC (Lower is Better)")
    p_icc  <- style_metric_plot(metrics_df, "ICC", "ICC (Intraclass Correlation)", is_icc = TRUE)
    p_pcv  <- style_metric_plot(pcv_df, "PCV", "PCV (Proportional Change in Variance)", is_pcv = TRUE)

    # 3. Combine into a professional 2x2 grid (AIC, BIC, ICC, VPR)
    cowplot::plot_grid(
      p_aic,
      style_metric_plot(metrics_df, "BIC", "BIC (Lower is Better)"),
      p_icc,
      p_pcv,
      ncol = 2,
      labels = "AUTO",
      label_size = 14
    )
  })
  output$freq_map <- renderPlot({
    req(freq_res())
    res <- freq_res()

    # Extract Random Effects (BLUPs)
    re <- as.data.frame(ranef(res$m3)$cond$cluster_idx)
    re$ID <- rownames(re)
    colnames(re)[1] <- "intercept_re"

    df_clean <- processed_df()
    cluster_region_map <- df_clean %>%
      group_by(!!sym(input$v_cluster)) %>%
      summarise(Region_Name = first(!!sym(input$v_reg))) %>%
      mutate(ID = as.character(!!sym(input$v_cluster)))

    re_mapped <- re %>%
      left_join(cluster_region_map, by = "ID") %>%
      group_by(Region_Name) %>%
      summarise(mean_random_effect = mean(intercept_re, na.rm = TRUE)) %>%
      mutate(Match_Key = str_to_upper(str_trim(as.character(Region_Name))))

    shp <- adal_spatial()
    map_df <- shp %>% left_join(re_mapped, by = "Match_Key")

    ggplot(map_df) +
      geom_sf(aes(fill = mean_random_effect), color = "white", size = 0.2) +
      scale_fill_gradient2(low = "blue", mid = "white", high = "red", name = "Random Effect") +
      theme_void() +
      labs(title = "Spatial Random Effects Map",
           subtitle = "Unexplained cluster-level variation aggregated to administrative level")
  })

  output$cov_corr_plot <- renderPlot({
    vars <- c(input$v_indiv_cont, input$v_indiv_fact,
              input$v_comm_cont, input$v_comm_fact)

    validate(
      need(length(vars) >= 2, "Please select at least 2 covariates in Setup.")
    )

    df_corr <- processed_df() %>%
      dplyr::select(all_of(vars)) %>% # Added dplyr:: prefix
      mutate(across(everything(), ~as.numeric(as.factor(.)))) %>%
      na.omit()

    M <- cor(df_corr)

    corrplot(M, method = "circle", type = "upper",
             order = "hclust", tl.col = "black", tl.srt = 45,
             addCoef.col = "black", number.cex = 0.7,
             col = brewer.pal(n = 8, name = "RdYlBu"),
             title = "Covariate Correlation Matrix",
             mar = c(0,0,2,0))
  })
  ###############################################################################################################
  # Dynamic UI for selecting waves based on the time variable
  output$wave_selector_1 <- renderUI({
    req(input$v_time, processed_df())
    choices <- unique(processed_df()[[input$v_time]])
    selectInput("wave1", "Reference Group (Wave 1):", choices = choices)
  })

  output$wave_selector_2 <- renderUI({
    req(input$v_time, processed_df())
    choices <- unique(processed_df()[[input$v_time]])
    selectInput("wave2", "Comparison Group (Wave 2):", choices = choices, selected = choices[length(choices)])
  })

  # Decomposition Engine
  # --- Tab 13: Spatiotemporal Decomposition Logic (Upgraded) ---
  decomp_res <- eventReactive(input$run_decomp, {
    req(input$wave1, input$wave2, input$v_dv)

    withProgress(message = 'Calculating Spatiotemporal Decomposition...', value = 0, {
      df <- processed_df()
      all_preds <- unique(c(input$v_indiv_cont, input$v_comm_cont, input$v_indiv_fact, input$v_comm_fact))

      incProgress(0.2, detail = "Preparing Waves")
      g1_data <- df %>% filter(!!sym(input$v_time) == input$wave1) %>% drop_na(dv_final, all_of(all_preds))
      g2_data <- df %>% filter(!!sym(input$v_time) == input$wave2) %>% drop_na(dv_final, all_of(all_preds))

      # Table 3: Rates by Category
      calc_st_t3 <- function(data) {
        res_list <- list()
        for(p in all_preds) {
          res_list[[p]] <- data %>% group_by(Category = as.character(haven::as_factor(!!sym(p)))) %>%
            summarise(Rate = weighted.mean(dv_final, !!sym(input$v_wgt), na.rm=TRUE), .groups='drop') %>%
            mutate(Variable = p)
        }
        do.call(rbind, res_list)
      }
      t3_final <- full_join(calc_st_t3(g1_data), calc_st_t3(g2_data), by=c("Variable", "Category"),
                            suffix=c(paste0("_", input$wave1), paste0("_", input$wave2))) %>%
        filter(!is.na(Category))

      incProgress(0.4, detail = "Fitting Models")
      valid_preds <- all_preds[sapply(all_preds, function(p) length(unique(g1_data[[p]])) > 1 && length(unique(g2_data[[p]])) > 1)]
      form <- as.formula(paste("dv_final ~", paste(valid_preds, collapse = " + ")))
      m1 <- glm(form, data = g1_data, family = binomial(link = "logit"))
      m2 <- glm(form, data = g2_data, family = binomial(link = "logit"))

      # Math for Inference
      b1 <- coef(m1); v1 <- vcov(m1); X1 <- model.matrix(m1); means1 <- colMeans(X1)
      b2 <- coef(m2); v2 <- vcov(m2); X2 <- model.matrix(m2); means2 <- colMeans(X2)
      diff_means <- means2 - means1

      endow_vec <- diff_means * b1
      coeff_vec <- means1 * (b2 - b1)

      y1 <- mean(g1_data$dv_final); y2 <- mean(g2_data$dv_final); gap <- y2 - y1
      se_e <- sqrt(drop(t(as.matrix(diff_means)) %*% v1 %*% as.matrix(diff_means)))
      se_c <- sqrt(drop(t(as.matrix(means1)) %*% (v1 + v2) %*% as.matrix(means1)))

      # Summary (Table 4)
      summary_df <- data.frame(
        Component = c(paste("Wave", input$wave1), paste("Wave", input$wave2), "Total Change", "Endowment (Explained)", "Coefficient (Unexplained)"),
        Value = c(y1, y2, gap, sum(endow_vec), sum(coeff_vec)),
        CI_Lower = c(NA, NA, NA, sum(endow_vec) - 1.96*se_e, sum(coeff_vec) - 1.96*se_c),
        CI_Upper = c(NA, NA, NA, sum(endow_vec) + 1.96*se_e, sum(coeff_vec) + 1.96*se_c),
        P_value = c(NA, NA, NA, 2*(1-pnorm(abs(sum(endow_vec)/se_e))), 2*(1-pnorm(abs(sum(coeff_vec)/se_c))))
      )

      # Detailed (Table 5/6)
      # 6. Detailed (Table 5/6) - AGGREGATED VERSION
      raw_detailed <- data.frame(Variable = names(b1)) %>%
        filter(Variable != "(Intercept)") %>%
        mutate(
          Endowment = as.numeric(endow_vec[match(Variable, names(b1))]),
          E_SE = sqrt(diff_means[match(Variable, names(b1))]^2 * diag(v1)[match(Variable, names(b1))]),
          Coefficient = as.numeric(coeff_vec[match(Variable, names(b1))]),
          C_SE = sqrt(means1[match(Variable, names(b1))]^2 * (diag(v1)[match(Variable, names(b1))] + diag(v2)[match(Variable, names(b1))])),
          Pct_Share = (Endowment / gap) * 100
        )

      # Map dummy variables back to parent variables
      raw_detailed$ParentVar <- NA
      for(p in all_preds) {
        raw_detailed$ParentVar[grepl(paste0("^", p), raw_detailed$Variable)] <- p
      }
      # Handle any variables that didn't match (fallback)
      raw_detailed$ParentVar[is.na(raw_detailed$ParentVar)] <- raw_detailed$Variable[is.na(raw_detailed$ParentVar)]

      detailed_df <- raw_detailed %>%
        group_by(Variable = ParentVar) %>%
        summarise(
          Endowment = sum(Endowment, na.rm=TRUE),
          E_SE = sqrt(sum(E_SE^2, na.rm=TRUE)),
          Coefficient = sum(Coefficient, na.rm=TRUE),
          C_SE = sqrt(sum(C_SE^2, na.rm=TRUE)),
          Pct_Share = sum(Pct_Share, na.rm=TRUE),
          .groups = 'drop'
        ) %>%
        mutate(
          E_Pval = 2*(1-pnorm(abs(Endowment/E_SE))),
          E_CI = paste0("(", round(Endowment - 1.96*E_SE, 3), ", ", round(Endowment + 1.96*E_SE, 3), ")"),
          C_Pval = 2*(1-pnorm(abs(Coefficient/C_SE))),
          C_CI = paste0("(", round(Coefficient - 1.96*C_SE, 3), ", ", round(Coefficient + 1.96*C_SE, 3), ")")
        )

      plot_df <- bind_rows(
        detailed_df %>% dplyr::select(Variable, Est=Endowment, SE=E_SE, P=E_Pval) %>% mutate(Type="Endowment (Explained)"),
        detailed_df %>% dplyr::select(Variable, Est=Coefficient, SE=C_SE, P=C_Pval) %>% mutate(Type="Coefficient (Unexplained)")
      )

      incProgress(0.4, detail = "Finalizing")

      # UPDATE THIS LIST:
      list(
        summary = summary_df,
        detail = detailed_df,
        plot_data = plot_df,
        table3 = t3_final,
        # ADD THIS LINE HERE:
        table2 = calc_descriptive_table2(g1_data, g2_data, input$wave1, input$wave2, all_preds, input$v_wgt),
        wave1 = input$wave1,
        wave2 = input$wave2
      )
    })
  })

  # Render Outputs for Tab 13
  output$st_table3 <- renderTable({ req(decomp_res()); decomp_res()$table3 }, digits = 4)
  output$st_decomp_summary_table <- renderTable({ req(decomp_res()); decomp_res()$summary }, digits = 4, na = "-")
  output$st_decomp_detail_table <- renderTable({
    req(decomp_res())
    decomp_res()$detail %>% dplyr::select(Variable, Endowment, E_CI, E_Pval, Coefficient, C_CI, C_Pval)
  }, digits = 4)

  output$st_figure5_plot <- renderPlot({
    req(decomp_res())
    df <- decomp_res()$plot_data
    ggplot(df, aes(x = Est, y = reorder(Variable, Est))) +
      geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +
      geom_errorbarh(aes(xmin = Est - 1.96*SE, xmax = Est + 1.96*SE), height = 0.2, color = 'gray') +
      geom_point(aes(color = P < 0.05), size = 3) +
      facet_wrap(~ Type, scales = "free_x") +
      scale_color_manual(values = c('TRUE' = '#c0392b', 'FALSE' = '#2c3e50'), name = 'Significant') +
      theme_minimal(base_size = 14) + labs(x = 'Contribution to Change', y = '', title = "Figure 5: Spatiotemporal Decomposition")
  })



  ################################################################################################
  # --- Spatio-Temporal Interaction Logic ---
  st_int_res <- eventReactive(input$run_st_int, {
    req(processed_df(), input$v_reg, input$v_time)

    df <- processed_df() %>%
      drop_na(dv_final, !!sym(input$v_reg), !!sym(input$v_time))

    # Create clean factors for the interaction term
    df$region_idx <- as.factor(df[[input$v_reg]])
    df$year_idx <- as.factor(df[[input$v_time]])
    df$year_num <- as.numeric(as.character(df[[input$v_time]]))

    preds <- unique(c(input$v_indiv_cont, input$v_indiv_fact, input$v_comm_cont, input$v_comm_fact))

    # The Interaction Formula
    pred_part <- if(length(preds) > 0) paste0(paste(preds, collapse = " + "), " + ") else ""
    form <- as.formula(paste0("dv_final ~ ", pred_part,
                              "year_num + (1 | region_idx) + (1 | year_idx) + (1 | region_idx:year_idx)"))

    withProgress(message = 'Fitting Spatio-Temporal Interaction...', value = 0.5, {
      model <- glmmTMB(form, data = df, family = binomial(link="logit"))
    })

    # Extract Interaction Effects
    re <- as.data.frame(ranef(model)$cond$`region_idx:year_idx`)
    re$ID <- rownames(re)

    # Get the interaction for the most recent year
    latest_year <- max(as.character(df$year_idx))

    re_latest <- re %>%
      filter(grepl(latest_year, ID)) %>%
      mutate(Match_Key = str_to_upper(str_trim(str_replace(ID, paste0(":", latest_year), "")))) %>%
      rename(interaction_effect = `(Intercept)`)

    return(list(model = model, interaction_data = re_latest))
  })

  output$st_int_map <- renderPlot({
    req(st_int_res())
    re_data <- st_int_res()$interaction_data
    shp <- adal_spatial() %>% dplyr::select(Match_Key, geometry)
    map_df <- shp %>% left_join(re_data, by = "Match_Key")

    ggplot(map_df) +
      geom_sf(aes(fill = interaction_effect), color = "white", size = 0.1) +
      scale_fill_gradient2(low = "#2980b9", mid = "white", high = "#c0392b", name = "Interaction") +
      theme_void() +
      labs(title = "Emerging Hotspots Map", subtitle = "Deviation from the national temporal trend")
  })

  output$st_int_metrics <- renderTable({
    req(st_int_res())
    m <- st_int_res()$model
    vc <- tryCatch(VarCorr(m), error = function(e) NULL)
    data.frame(
      Metric = c("Observations", "AIC", "BIC", "Log-likelihood", "Fixed-effect p-values", "Region-time interaction"),
      Value = c(
        as.character(stats::nobs(m)),
        round(AIC(m), 3),
        round(BIC(m), 3),
        round(as.numeric(logLik(m)), 3),
        "Reported in the coefficient table",
        ifelse(is.null(vc), "Estimated", "Estimated")
      )
    )
  }, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m", width = "100%")

  output$st_int_coef_table <- renderTable({
    req(st_int_res())
    tidy_glmmTMB_fixed_apa(st_int_res()$model, exponentiate = TRUE)
  }, digits = 3, striped = TRUE, bordered = TRUE, hover = TRUE, spacing = "m", width = "100%")
  ###############################################################################################################
  # --- Descriptive Statistics Table Logic (Fixes Visibility & Strict Requirements) ---

  desc_table_data <- reactive({
    # 1. Basic Requirement: Data must be uploaded
    req(raw_data())

    # 2. Collect variables (Handle NULLs)
    cont_vars <- c(input$v_indiv_cont, input$v_comm_cont)
    cat_vars <- c(input$v_indiv_fact, input$v_comm_fact)

    # 3. Validation: If NO variables are selected at all, return a placeholder
    if (length(cont_vars) == 0 && length(cat_vars) == 0) {
      return(data.frame(Variable = "Please select variables", Category = "in Data Setup tab", Statistic = ""))
    }

    df <- raw_data()

    # 4. Weight Handling (Robust check)
    if (!is.null(input$v_wgt) && input$v_wgt != "") {
      if(input$v_wgt %in% names(df)) {
        wt_raw <- df[[input$v_wgt]]
        # Handle DHS scaling (divide by 1M) or standard weight
        wt <- if(input$is_dhs_wgt) wt_raw / 1000000 else wt_raw
      } else {
        wt <- rep(1, nrow(df))
      }
    } else {
      wt <- rep(1, nrow(df))
    }

    # Clean weights: replace NAs with 0 to prevent calculation errors
    wt[is.na(wt)] <- 0

    table_rows <- list()

    # --- A. Continuous Variables (Mean/SD) ---
    for(v in cont_vars) {
      if(v %in% names(df)) {
        val <- as.numeric(df[[v]]) # Force numeric
        # Filter valid data (non-NA value, non-NA weight, positive weight)
        valid_idx <- !is.na(val) & !is.na(wt) & wt > 0

        if(sum(valid_idx) > 0) {
          val_c <- val[valid_idx]
          wt_c <- wt[valid_idx]

          w_mean <- weighted.mean(val_c, wt_c)
          w_sd <- sqrt(cov.wt(data.frame(val_c), wt = wt_c)$cov[1,1])

          table_rows[[length(table_rows) + 1]] <- data.frame(
            Variable = v,
            Category = "Mean (SD)",
            Statistic = paste0(round(w_mean, 2), " (", round(w_sd, 2), ")")
          )
        }
      }
    }

    # --- B. Categorical Variables (Freq/%) ---
    for(v in cat_vars) {
      if(v %in% names(df)) {
        # Use haven::as_factor to get LABELS (e.g., "Urban") instead of numbers
        val_factor <- haven::as_factor(df[[v]])

        # Create temp dataframe and filter
        tmp <- data.frame(cat = val_factor, w = wt) %>%
          filter(!is.na(cat), !is.na(w), w > 0)

        if(nrow(tmp) > 0) {
          agg <- tmp %>%
            group_by(cat) %>%
            summarise(n_w = sum(w), .groups = 'drop') %>%
            mutate(pct = (n_w / sum(n_w)) * 100)

          # Add Header Row for the Variable
          table_rows[[length(table_rows) + 1]] <- data.frame(
            Variable = as.character(v), Category = "", Statistic = ""
          )

          # Add Rows for each Category
          for(i in 1:nrow(agg)) {
            table_rows[[length(table_rows) + 1]] <- data.frame(
              Variable = "",
              Category = as.character(agg$cat[i]),
              Statistic = paste0(round(agg$n_w[i], 0), " (", round(agg$pct[i], 1), "%)")
            )
          }
        }
      }
    }

    # Final check: if no rows generated (e.g. all data was NA)
    if(length(table_rows) == 0) return(data.frame(Variable = "No valid data found", Category="-", Statistic="-"))

    do.call(rbind, table_rows)
  })

  # 2. Render the Table (Updated for Visibility and Style)
  output$desc_stats_table <- renderTable({
    # Get the data
    tbl <- desc_table_data()

    # Safety check: if the result is just a message, return it
    if(is.null(tbl) || nrow(tbl) == 0) return(data.frame(Status = "Please select variables in Data Setup tab."))

    # Return the table
    tbl
  },
  # --- STYLING ARGUMENTS TO MATCH MULTILEVEL TABLE ---
  striped = TRUE,
  bordered = TRUE,  # Adds the borders
  hover = TRUE,
  width = "100%",
  align = 'l',
  spacing = 'm'     # Adjusts row height
  )


  ####################################################################################################################
  # 2. Core Engine
  # --- Tab 14: Regional Decomposition Logic (Paper Replication) ---

  # 1. UI Selectors for Groups
  output$risk_selector_1 <- renderUI({
    if(isTRUE(input$auto_lisa)) return(NULL)
    req(input$v_risk_group, processed_df())
    choices <- unique(as.character(processed_df()[[input$v_risk_group]]))
    selectInput("risk1", "Reference Group (G1):", choices = choices)
  })

  output$risk_selector_2 <- renderUI({
    if(isTRUE(input$auto_lisa)) return(NULL)
    req(input$v_risk_group, processed_df())
    choices <- unique(as.character(processed_df()[[input$v_risk_group]]))
    selectInput("risk2", "Comparison Group (G2):", choices = choices)
  })

  # 2. The Core Decomposition Engine
  spatial_decomp_res <- eventReactive(input$run_spatial_decomp, {
    req(input$v_dv)
    df <- processed_df()
    df <- df %>% mutate(Match_Key = str_to_upper(str_trim(as.character(!!sym(input$v_reg)))))
    all_preds <- unique(c(input$v_indiv_cont, input$v_comm_cont, input$v_indiv_fact, input$v_comm_fact))

    # 1. Group Selection
    if(input$auto_lisa) {
      shp <- adal_spatial()
      nb <- poly2nb(shp, queen = TRUE)
      if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(shp)), k = 1))
      lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
      loc_m <- localmoran(shp$outcome, lw, zero.policy = TRUE)
      m_avg <- mean(shp$outcome, na.rm=TRUE)
      shp$lisa_group <- case_when(loc_m[,5] <= 0.05 & shp$outcome >= m_avg ~ "Hot Spot",
                                  loc_m[,5] <= 0.05 & shp$outcome <= m_avg ~ "Cold Spot", TRUE ~ "Other")
      df <- df %>% left_join(st_drop_geometry(shp) %>% dplyr::select(Match_Key, lisa_group), by = "Match_Key")
      g1_label <- "Cold Spot"; g2_label <- "Hot Spot"
      g1_data <- df %>% filter(lisa_group == "Cold Spot"); g2_data <- df %>% filter(lisa_group == "Hot Spot")
    } else {
      g1_label <- input$risk1; g2_label <- input$risk2
      g1_data <- df %>% filter(!!sym(input$v_risk_group) == input$risk1); g2_data <- df %>% filter(!!sym(input$v_risk_group) == input$risk2)
    }

    # 2. Table 3: Rates by Category (Cleaned)
    calc_t3 <- function(data) {
      res_list <- list()
      for(p in all_preds) {
        res_list[[p]] <- data %>% group_by(Category = as.character(haven::as_factor(!!sym(p)))) %>%
          summarise(Rate = weighted.mean(dv_final, !!sym(input$v_wgt), na.rm=TRUE), .groups='drop') %>%
          mutate(Variable = p)
      }
      do.call(rbind, res_list)
    }
    t3_final <- full_join(calc_t3(g1_data), calc_t3(g2_data), by=c("Variable", "Category"),
                          suffix=c(paste0("_", g1_label), paste0("_", g2_label))) %>%
      filter(!is.na(Category))

    # 3. Models
    g1_m <- g1_data %>% drop_na(dv_final, all_of(all_preds))
    g2_m <- g2_data %>% drop_na(dv_final, all_of(all_preds))
    valid_preds <- all_preds[sapply(all_preds, function(p) length(unique(g1_m[[p]])) > 1 && length(unique(g2_m[[p]])) > 1)]
    form <- as.formula(paste("dv_final ~", paste(valid_preds, collapse = " + ")))
    m1 <- glm(form, data = g1_m, family = binomial(link = "logit"))
    m2 <- glm(form, data = g2_m, family = binomial(link = "logit"))

    # 4. Math
    b1 <- coef(m1); v1 <- vcov(m1); X1 <- model.matrix(m1); means1 <- colMeans(X1)
    b2 <- coef(m2); v2 <- vcov(m2); X2 <- model.matrix(m2); means2 <- colMeans(X2)
    diff_means <- means2 - means1

    endow_vec <- diff_means * b1
    coeff_vec <- means1 * (b2 - b1)

    y1 <- mean(g1_m$dv_final); y2 <- mean(g2_m$dv_final); gap <- y2 - y1
    se_e <- sqrt(drop(t(as.matrix(diff_means)) %*% v1 %*% as.matrix(diff_means)))
    se_c <- sqrt(drop(t(as.matrix(means1)) %*% (v1 + v2) %*% as.matrix(means1)))

    # 5. Summary (Table 4)
    summary_df <- data.frame(
      Component = c(paste(g1_label, "Mean"), paste(g2_label, "Mean"), "Total Gap", "Endowment (Explained)", "Coefficient (Unexplained)"),
      Value = c(y1, y2, gap, sum(endow_vec), sum(coeff_vec)),
      CI_Lower = c(NA, NA, NA, sum(endow_vec) - 1.96*se_e, sum(coeff_vec) - 1.96*se_c),
      CI_Upper = c(NA, NA, NA, sum(endow_vec) + 1.96*se_e, sum(coeff_vec) + 1.96*se_c),
      P_value = c(NA, NA, NA, 2*(1-pnorm(abs(sum(endow_vec)/se_e))), 2*(1-pnorm(abs(sum(coeff_vec)/se_c))))
    )

    # 6. Detailed (Table 5/6)
    # 6. Detailed (Table 5/6) - AGGREGATED VERSION
    raw_detailed_sp <- data.frame(Variable = names(b1)) %>%
      filter(Variable != "(Intercept)") %>%
      mutate(
        Endowment = as.numeric(endow_vec[match(Variable, names(b1))]),
        E_SE = sqrt(diff_means[match(Variable, names(b1))]^2 * diag(v1)[match(Variable, names(b1))]),
        Coefficient = as.numeric(coeff_vec[match(Variable, names(b1))]),
        C_SE = sqrt(means1[match(Variable, names(b1))]^2 * (diag(v1)[match(Variable, names(b1))] + diag(v2)[match(Variable, names(b1))])),
        Pct_Share = (Endowment / gap) * 100
      )

    raw_detailed_sp$ParentVar <- NA
    for(p in all_preds) {
      raw_detailed_sp$ParentVar[grepl(paste0("^", p), raw_detailed_sp$Variable)] <- p
    }
    raw_detailed_sp$ParentVar[is.na(raw_detailed_sp$ParentVar)] <- raw_detailed_sp$Variable[is.na(raw_detailed_sp$ParentVar)]

    detailed_df <- raw_detailed_sp %>%
      group_by(Variable = ParentVar) %>%
      summarise(
        Endowment = sum(Endowment, na.rm=TRUE),
        E_SE = sqrt(sum(E_SE^2, na.rm=TRUE)),
        Coefficient = sum(Coefficient, na.rm=TRUE),
        C_SE = sqrt(sum(C_SE^2, na.rm=TRUE)),
        Pct_Share = sum(Pct_Share, na.rm=TRUE),
        .groups = 'drop'
      ) %>%
      mutate(
        E_Pval = 2*(1-pnorm(abs(Endowment/E_SE))),
        E_CI = paste0("(", round(Endowment - 1.96*E_SE, 3), ", ", round(Endowment + 1.96*E_SE, 3), ")"),
        C_Pval = 2*(1-pnorm(abs(Coefficient/C_SE))),
        C_CI = paste0("(", round(Coefficient - 1.96*C_SE, 3), ", ", round(Coefficient + 1.96*C_SE, 3), ")")
      )

    plot_df <- bind_rows(
      detailed_df %>% dplyr::select(Variable, Est=Endowment, SE=E_SE, P=E_Pval) %>% mutate(Type="Endowment (Explained)"),
      detailed_df %>% dplyr::select(Variable, Est=Coefficient, SE=C_SE, P=C_Pval) %>% mutate(Type="Coefficient (Unexplained)")
    )

    list(
      summary = summary_df,
      detail = detailed_df,
      plot_data = plot_df,
      table3 = t3_final,
      # ADD THIS LINE HERE:
      table2 = calc_descriptive_table2(g1_data, g2_data, g1_label, g2_label, all_preds, input$v_wgt),
      g1_name = g1_label,
      g2_name = g2_label
    )
  })

  # 3. Render Outputs (Assigning to UI IDs)
  output$spatial_table3 <- renderTable({ req(spatial_decomp_res()); spatial_decomp_res()$table3 }, digits = 4)

  output$spatial_decomp_summary_table <- renderTable({ req(spatial_decomp_res()); spatial_decomp_res()$summary }, digits = 4, na = "-")

  output$spatial_decomp_detail_table <- renderTable({
    req(spatial_decomp_res())
    spatial_decomp_res()$detail %>% dplyr::select(Variable, Endowment, E_CI, E_Pval, Coefficient, C_CI, C_Pval)
  }, digits = 4)

  output$spatial_figure5_plot <- renderPlot({
    req(spatial_decomp_res())
    df <- spatial_decomp_res()$plot_data

    ggplot(df, aes(x = Est, y = reorder(Variable, Est))) +
      geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +
      geom_errorbarh(aes(xmin = Est - 1.96*SE, xmax = Est + 1.96*SE), height = 0.2, color = 'gray') +
      geom_point(aes(color = P < 0.05), size = 3) +
      facet_wrap(~ Type, scales = "free_x") +
      scale_color_manual(values = c('TRUE' = '#c0392b', 'FALSE' = '#2c3e50'), name = 'Significant') +
      theme_minimal(base_size = 14) +
      labs(x = 'Contribution to Gap', y = '', title = "Figure 5: Decomposition of Regional Gap")
  })



  # --- Render Table 2 ---
  output$st_table2 <- renderTable({ req(decomp_res()); decomp_res()$table2 }, digits = 1)
  output$spatial_table2 <- renderTable({ req(spatial_decomp_res()); spatial_decomp_res()$table2 }, digits = 1)

  # --- Render Percentage Contribution Plots ---
  output$st_pct_plot <- renderPlot({
    req(decomp_res())
    df <- decomp_res()$summary[4:5,]
    df$Pct <- (df$Value / sum(df$Value)) * 100
    ggplot(df, aes(x = Component, y = Pct, fill = Component)) +
      geom_bar(stat = "identity", width = 0.6) +
      geom_text(aes(label = paste0(round(Pct, 1), "%")), vjust = -0.5, fontface="bold") +
      scale_fill_manual(values = c("#2c3e50", "#e74c3c")) +
      theme_minimal() + labs(title = "Total Contribution (%)", y = "Percentage Share", x = "")
  })

  output$spatial_pct_plot <- renderPlot({
    req(spatial_decomp_res())
    df <- spatial_decomp_res()$summary[4:5,]
    df$Pct <- (df$Value / sum(df$Value)) * 100
    ggplot(df, aes(x = Component, y = Pct, fill = Component)) +
      geom_bar(stat = "identity", width = 0.6) +
      geom_text(aes(label = paste0(round(Pct, 1), "%")), vjust = -0.5, fontface="bold") +
      scale_fill_manual(values = c("#2c3e50", "#e74c3c")) +
      theme_minimal() + labs(title = "Total Contribution (%)", y = "Percentage Share", x = "")
  })


  ########################################################################################################
  # --- 1. Joint Model Engine (glmmTMB) ---
  joint_res <- eventReactive(input$run_joint, {
    req(input$joint_dv1, input$joint_dv2, processed_df())
    df <- processed_df()
    preds <- unique(c(input$v_indiv_cont, input$v_indiv_fact, input$v_comm_cont, input$v_comm_fact))

    # Reshape to Long
    df_long <- df %>%
      dplyr::select(all_of(preds),
                    outcome1 = !!sym(input$joint_dv1),
                    outcome2 = !!sym(input$joint_dv2),
                    cluster_idx = !!sym(input$v_cluster)) %>%
      pivot_longer(cols = c(outcome1, outcome2), names_to = "outcome_type", values_to = "joint_value") %>%
      mutate(outcome_type = as.factor(outcome_type), cluster_idx = as.factor(cluster_idx))

    form <- as.formula(paste("joint_value ~ 0 + outcome_type + outcome_type:(",
                             paste(preds, collapse = " + "), ") + (0 + outcome_type | cluster_idx)"))

    withProgress(message = 'Fitting Joint Model...', value = 0.5, {
      glmmTMB(form, data = df_long, family = binomial(link = "logit"))
    })
  })

  # --- 2. Joint Spatial Aggregator (For Moran's I) ---
  # This calculates the % of people in each region who have BOTH conditions
  joint_spatial_data <- reactive({
    req(input$run_joint)
    df <- processed_df()
    df$wgt_f <- if(input$is_dhs_wgt) df[[input$v_wgt]]/1000000 else as.numeric(df[[input$v_wgt]])

    # Create "Both" indicator
    df$both_conditions <- ifelse(df[[input$joint_dv1]] == 1 & df[[input$joint_dv2]] == 1, 1, 0)

    stats_joint <- df %>%
      group_by(Match_Key = str_to_upper(str_trim(as.character(!!sym(input$v_reg))))) %>%
      summarise(joint_rate = weighted.mean(both_conditions, wgt_f, na.rm=TRUE))

    adal_spatial() %>% left_join(stats_joint, by = "Match_Key")
  })

  # --- 3. Shared Variance Correlation Box ---
  output$joint_corr_box <- renderInfoBox({
    req(joint_res())
    vc <- VarCorr(joint_res())
    corr_val <- attr(vc$cond$cluster_idx, "correlation")[1,2]
    infoBox("Shared Correlation", round(corr_val, 3), icon = icon("link"), color = "purple", fill = TRUE)
  })

  # --- 4. Global Moran's I for Joint Outcome ---
  output$joint_moran_box <- renderInfoBox({
    df <- joint_spatial_data()
    nb <- poly2nb(df, queen = TRUE)
    if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

    m_test <- moran.test(df$joint_rate, lw, zero.policy = TRUE)
    infoBox("Joint Moran's I", round(m_test$estimate[1], 3), icon = icon("th-large"), color = "maroon", fill = TRUE)
  })

  output$joint_moran_text <- renderPrint({
    df <- joint_spatial_data()
    nb <- poly2nb(df, queen = TRUE)
    if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
    moran.test(df$joint_rate, lw, zero.policy = TRUE)
  })

  # --- 5. FIXED Forest Plot ---
  output$joint_forest_plot <- renderPlot({
    req(joint_res())
    # Tidy the model and clean up the term names for plotting
    res_df <- broom.mixed::tidy(joint_res(), effects = "fixed", conf.int = TRUE) %>%
      filter(!term %in% c("outcome_typeoutcome1", "outcome_typeoutcome2")) %>%
      mutate(
        # Clean names: "outcome_typeoutcome1:age" -> "age (Outcome 1)"
        Outcome = ifelse(str_detect(term, "outcome1"), "Outcome 1", "Outcome 2"),
        Variable = str_replace(term, "outcome_typeoutcome[12]:", "")
      )

    ggplot(res_df, aes(x = estimate, y = Variable, color = Outcome)) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
      geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.3, position = position_dodge(width = 0.5)) +
      geom_point(size = 3, position = position_dodge(width = 0.5)) +
      scale_color_manual(values = c("Outcome 1" = "#2c3e50", "Outcome 2" = "#e67e22")) +
      theme_minimal(base_size = 14) +
      labs(title = "Joint Determinants (Log-Odds)", x = "Effect Size", y = "")
  })

  # --- 6. LISA Map for Joint Outcome ---
  output$joint_lisa_map <- renderPlot({
    df <- joint_spatial_data()
    nb <- poly2nb(df, queen = TRUE)
    if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

    loc_m <- localmoran(df$joint_rate, lw, zero.policy = TRUE)
    df$Pr <- loc_m[,5]
    df$lag_joint <- lag.listw(lw, df$joint_rate, zero.policy = TRUE)
    m_avg <- mean(df$joint_rate, na.rm=TRUE)

    df$lisa <- case_when(
      df$Pr <= 0.05 & df$joint_rate >= m_avg & df$lag_joint >= m_avg ~ "High-High (Hotspot)",
      df$Pr <= 0.05 & df$joint_rate <= m_avg & df$lag_joint <= m_avg ~ "Low-Low (Coldspot)",
      TRUE ~ "Not Significant"
    )

    ggplot(df) + geom_sf(aes(fill = lisa), color = "black") +
      scale_fill_manual(values = c("High-High (Hotspot)"="red", "Low-Low (Coldspot)"="blue", "Not Significant"="lightgray")) +
      theme_void() + labs(title = "Joint LISA Cluster Map", fill = "Cluster Type")
  })

  # --- 7. Table ---
  output$joint_coef_table <- renderTable({
    req(joint_res())
    res <- broom.mixed::tidy(joint_res(), effects = "fixed", conf.int = TRUE)
    # Convert to AOR
    res %>% mutate(across(c(estimate, conf.low, conf.high), exp))
  }, digits = 3)

  ##############################################################################################################
  # 1. Status Messenger
  jd_status <- reactiveVal("Waiting for inputs...")

  # 2. Selectors
  output$jd_risk_selector_1 <- renderUI({
    if(is.null(input$jd_group) || input$jd_group == "") return(NULL)
    choices <- unique(as.character(processed_df()[[input$jd_group]]))
    selectInput("jd_risk1", "Reference Group (G1):", choices = choices)
  })

  output$jd_risk_selector_2 <- renderUI({
    if(is.null(input$jd_group) || input$jd_group == "") return(NULL)
    choices <- unique(as.character(processed_df()[[input$jd_group]]))
    selectInput("jd_risk2", "Comparison Group (G2):", choices = choices, selected = choices[length(choices)])
  })

  # 3. The Engine
  # --- 16. JOINT REGIONAL DECOMPOSITION ENGINE (UPGRADED) ---

  # Engine
  jd_res <- eventReactive(input$run_joint_decomp, {
    validate(
      need(input$jd_dv1 != "" && input$jd_dv2 != "", "Select both outcomes"),
      need(input$jd_group != "", "Select grouping variable"),
      need(input$jd_risk1 != "" && input$jd_risk2 != "", "Select groups"),
      need(input$jd_risk1 != input$jd_risk2, "Groups must be different")
    )

    withProgress(message = 'Calculating Joint Decomposition...', value = 0, {
      df <- processed_df()

      # Robust Outcome Prep
      prep_y <- function(vec) {
        v <- as.numeric(as.character(vec))
        if (all(v %in% c(1, 2), na.rm = TRUE)) v <- v - 1
        if (max(v, na.rm=T) > 1) v <- ifelse(v == max(v, na.rm=T), 1, 0)
        return(v)
      }
      df$joint_y <- ifelse(prep_y(df[[input$jd_dv1]]) == 1 & prep_y(df[[input$jd_dv2]]) == 1, 1, 0)

      all_preds <- unique(c(input$v_indiv_cont, input$v_comm_cont, input$v_indiv_fact, input$v_comm_fact))
      all_preds <- setdiff(all_preds, input$jd_group)

      g1_data <- df %>% filter(as.character(!!sym(input$jd_group)) == as.character(input$jd_risk1)) %>% drop_na(joint_y)
      g2_data <- df %>% filter(as.character(!!sym(input$jd_group)) == as.character(input$jd_risk2)) %>% drop_na(joint_y)

      # Table 3: Rates by Category
      calc_jd_t3 <- function(data) {
        res_list <- list()
        for(p in all_preds) {
          res_list[[p]] <- data %>% group_by(Category = as.character(haven::as_factor(!!sym(p)))) %>%
            summarise(Rate = weighted.mean(joint_y, !!sym(input$v_wgt), na.rm=TRUE), .groups='drop') %>%
            mutate(Variable = p)
        }
        do.call(rbind, res_list)
      }
      t3_final <- full_join(calc_jd_t3(g1_data), calc_jd_t3(g2_data), by=c("Variable", "Category"),
                            suffix=c(paste0("_", input$jd_risk1), paste0("_", input$jd_risk2))) %>%
        filter(!is.na(Category))

      # Models
      valid_preds <- all_preds[sapply(all_preds, function(p) length(unique(na.omit(g1_data[[p]]))) > 1 && length(unique(na.omit(g2_data[[p]]))) > 1)]
      form <- as.formula(paste("joint_y ~", paste(valid_preds, collapse = " + ")))
      m1 <- glm(form, data = g1_data, family = binomial(link = "logit"))
      m2 <- glm(form, data = g2_data, family = binomial(link = "logit"))

      # Math & Alignment
      all_terms <- union(names(coef(m1)), names(coef(m2)))
      b1 <- setNames(numeric(length(all_terms)), all_terms); b1[names(coef(m1))] <- coef(m1)
      b2 <- setNames(numeric(length(all_terms)), all_terms); b2[names(coef(m2))] <- coef(m2)
      v1 <- matrix(0, length(all_terms), length(all_terms), dimnames=list(all_terms, all_terms))
      v1[rownames(vcov(m1)), colnames(vcov(m1))] <- as.matrix(vcov(m1))
      v2 <- matrix(0, length(all_terms), length(all_terms), dimnames=list(all_terms, all_terms))
      v2[rownames(vcov(m2)), colnames(vcov(m2))] <- as.matrix(vcov(m2))
      means1 <- setNames(numeric(length(all_terms)), all_terms); means1[colnames(model.matrix(m1))] <- colMeans(model.matrix(m1))
      means2 <- setNames(numeric(length(all_terms)), all_terms); means2[colnames(model.matrix(m2))] <- colMeans(model.matrix(m2))

      diff_means <- means2 - means1
      endow_vec <- diff_means * b1
      coeff_vec <- means1 * (b2 - b1)
      gap <- mean(g2_data$joint_y) - mean(g1_data$joint_y)

      # Detailed Aggregation (Summing dummies to Parent Variables)
      raw_det <- data.frame(Term = all_terms) %>% filter(Term != "(Intercept)") %>%
        mutate(Endowment = endow_vec[match(Term, all_terms)], Coefficient = coeff_vec[match(Term, all_terms)])

      raw_det$ParentVar <- NA
      for(p in all_preds) { raw_det$ParentVar[grepl(paste0("^", p), raw_det$Term)] <- p }

      detailed_agg <- raw_det %>% group_by(Variable = ParentVar) %>%
        summarise(Endowment = sum(Endowment), Coefficient = sum(Coefficient), .groups = 'drop') %>%
        mutate(Pct_Share = (Endowment / gap) * 100)

      list(
        summary = data.frame(Component = c("G1 Mean", "G2 Mean", "Total Gap", "Explained (Endowment)", "Unexplained (Coefficient)"),
                             Value = c(mean(g1_data$joint_y), mean(g2_data$joint_y), gap, sum(endow_vec), sum(coeff_vec))),
        table2 = calc_descriptive_table2(g1_data, g2_data, input$jd_risk1, input$jd_risk2, all_preds, input$v_wgt),
        table3 = t3_final,
        detail = detailed_agg
      )
    })
  })

  # Outputs for UI
  output$jd_summary_table <- renderTable({ jd_res()$summary }, digits = 4)
  output$jd_detail_table <- renderTable({ jd_res()$detail }, digits = 4)
  output$jd_forest_plot <- renderPlot({
    req(jd_res())
    ggplot(jd_res()$detail, aes(x = Endowment, y = reorder(Variable, Endowment))) +
      geom_vline(xintercept = 0, linetype = "dashed") +
      geom_point(size = 4, color = "#c0392b") + theme_minimal() + labs(title="Drivers of Joint Gap (Explained Part)", x="Contribution", y="")
  })
  output$jd_pct_plot <- renderPlot({
    req(jd_res())
    df <- jd_res()$summary[4:5,]
    df$Pct <- (df$Value / sum(df$Value)) * 100
    ggplot(df, aes(x = Component, y = Pct, fill = Component)) +
      geom_bar(stat = "identity", width = 0.6) + scale_fill_manual(values = c("#2c3e50", "#e74c3c")) +
      theme_minimal() + labs(title = "% Contribution to Joint Gap", y = "Percentage Share", x = "")
  })
  ##############################################################################################################
  # 1. Selectors for Tab 17
  output$jsd_wave_selector_1 <- renderUI({
    req(input$v_time, processed_df())
    choices <- unique(as.character(processed_df()[[input$v_time]]))
    selectInput("jsd_wave1", "Baseline Wave (T1):", choices = choices)
  })

  output$jsd_wave_selector_2 <- renderUI({
    req(input$v_time, processed_df())
    choices <- unique(as.character(processed_df()[[input$v_time]]))
    selectInput("jsd_wave2", "Endline Wave (T2):", choices = choices, selected = choices[length(choices)])
  })

  # 2. The Engine (JSD)
  jsd_res <- eventReactive(input$run_jsd, {
    validate(
      need(input$jsd_dv1 != "" && input$jsd_dv2 != "", "Select outcomes"),
      need(input$jsd_wave1 != input$jsd_wave2, "Waves must be different")
    )

    withProgress(message = 'Calculating Joint Spatiotemporal Decomp...', value = 0, {
      df <- processed_df()

      # Robust Outcome Prep
      prep_y <- function(vec) {
        v <- as.numeric(as.character(vec))
        if (all(v %in% c(1, 2), na.rm = TRUE)) v <- v - 1
        if (max(v, na.rm=T) > 1) v <- ifelse(v == max(v, na.rm=T), 1, 0)
        return(v)
      }
      df$joint_y <- ifelse(prep_y(df[[input$jsd_dv1]]) == 1 & prep_y(df[[input$jsd_dv2]]) == 1, 1, 0)

      all_preds <- unique(c(input$v_indiv_cont, input$v_comm_cont, input$v_indiv_fact, input$v_comm_fact))

      g1_data <- df %>% filter(as.character(!!sym(input$v_time)) == as.character(input$jsd_wave1)) %>% drop_na(joint_y)
      g2_data <- df %>% filter(as.character(!!sym(input$v_time)) == as.character(input$jsd_wave2)) %>% drop_na(joint_y)

      # Models & Robust Alignment
      valid_preds <- all_preds[sapply(all_preds, function(p) length(unique(na.omit(g1_data[[p]]))) > 1 && length(unique(na.omit(g2_data[[p]]))) > 1)]
      form <- as.formula(paste("joint_y ~", paste(valid_preds, collapse = " + ")))
      m1 <- glm(form, data = g1_data, family = binomial(link = "logit"))
      m2 <- glm(form, data = g2_data, family = binomial(link = "logit"))

      all_terms <- union(names(coef(m1)), names(coef(m2)))
      b1 <- setNames(numeric(length(all_terms)), all_terms); b1[names(coef(m1))] <- coef(m1)
      b2 <- setNames(numeric(length(all_terms)), all_terms); b2[names(coef(m2))] <- coef(m2)
      means1 <- setNames(numeric(length(all_terms)), all_terms); means1[colnames(model.matrix(m1))] <- colMeans(model.matrix(m1))
      means2 <- setNames(numeric(length(all_terms)), all_terms); means2[colnames(model.matrix(m2))] <- colMeans(model.matrix(m2))

      # Math
      endow_vec <- (means2 - means1) * b1
      coeff_vec <- means1 * (b2 - b1)
      gap <- mean(g2_data$joint_y) - mean(g1_data$joint_y)

      # Aggregation
      raw_det <- data.frame(Term = all_terms) %>% filter(Term != "(Intercept)") %>%
        mutate(Endowment = endow_vec[match(Term, all_terms)], Coefficient = coeff_vec[match(Term, all_terms)])
      raw_det$ParentVar <- NA
      for(p in all_preds) { raw_det$ParentVar[grepl(paste0("^", p), raw_det$Term)] <- p }

      detailed_agg <- raw_det %>% group_by(Variable = ParentVar) %>%
        summarise(Endowment = sum(Endowment), Coefficient = sum(Coefficient), .groups = 'drop') %>%
        mutate(Pct_Share = (Endowment / gap) * 100)

      list(
        summary = data.frame(Component = c("Wave 1 Mean", "Wave 2 Mean", "Total Change", "Explained (Endowment)", "Unexplained (Coefficient)"),
                             Value = c(mean(g1_data$joint_y), mean(g2_data$joint_y), gap, sum(endow_vec), sum(coeff_vec))),
        table2 = calc_descriptive_table2(g1_data, g2_data, input$jsd_wave1, input$jsd_wave2, all_preds, input$v_wgt),
        detail = detailed_agg
      )
    })
  })

  # Outputs
  output$jsd_summary_table <- renderTable({ jsd_res()$summary })
  output$jsd_table2 <- renderTable({ jsd_res()$table2 }, digits = 1)
  output$jsd_detail_table <- renderTable({ jsd_res()$detail }, digits = 4)
  output$jsd_forest_plot <- renderPlot({
    req(jsd_res())
    ggplot(jsd_res()$detail, aes(x = Endowment, y = reorder(Variable, Endowment))) +
      geom_vline(xintercept = 0, linetype = "dashed") +
      geom_point(size = 4, color = "#27ae60") + theme_minimal() + labs(title="Drivers of Joint Change Over Time")
  })


  ########################################################################################################
  # --- PAPERS 11-15: ADDITIONAL SPATIAL PUBLICATION MODULES ---

  service_points_raw <- reactive({
    req(input$p11_service_csv)
    readr::read_csv(input$p11_service_csv$datapath, show_col_types = FALSE)
  })

  output$p11_lat_ui <- renderUI({
    if (is.null(input$p11_service_csv)) return(helpText("No service-point file uploaded. The app will use the selected bivariate X/access variable as a proxy."))
    selectInput("p11_lat", "Latitude column:", choices = names(service_points_raw()))
  })

  output$p11_lon_ui <- renderUI({
    if (is.null(input$p11_service_csv)) return(NULL)
    selectInput("p11_lon", "Longitude column:", choices = names(service_points_raw()))
  })

  paper11_res <- eventReactive(input$run_p11, {
    req(adal_spatial())
    sf_df <- adal_spatial() %>% dplyr::filter(!is.na(outcome))
    cent <- sf::st_centroid(sf_df)
    method <- "Bivariate/access proxy"

    if (!is.null(input$p11_service_csv) && !is.null(input$p11_lat) && !is.null(input$p11_lon) &&
        input$p11_lat %in% names(service_points_raw()) && input$p11_lon %in% names(service_points_raw())) {
      pts_raw <- service_points_raw() %>%
        dplyr::mutate(.lat = suppressWarnings(as.numeric(.data[[input$p11_lat]])),
                      .lon = suppressWarnings(as.numeric(.data[[input$p11_lon]]))) %>%
        dplyr::filter(!is.na(.lat), !is.na(.lon))
      validate(need(nrow(pts_raw) > 0, "Service-points CSV has no valid latitude/longitude rows."))
      pts <- sf::st_as_sf(pts_raw, coords = c(".lon", ".lat"), crs = 4326, remove = FALSE) %>% sf::st_transform(sf::st_crs(sf_df))
      dmat <- sf::st_distance(cent, pts)
      sf_df$Access_Metric <- as.numeric(apply(dmat, 1, min)) / 1000
      sf_df$Access_Label <- "Distance to nearest service (km)"
      method <- "Nearest service distance"
    } else {
      validate(need("bi_var" %in% names(sf_df), "Upload a service-points CSV or select a bivariate X variable in Data Setup."))
      sf_df$Access_Metric <- as.numeric(sf_df$bi_var)
      sf_df$Access_Label <- paste0("Access proxy: ", input$v_wealth)
    }

    if (input$p11_desert_rule == "median") {
      cutoff <- stats::median(sf_df$Access_Metric, na.rm = TRUE)
      if (method == "Nearest service distance") {
        sf_df$Desert_Status <- ifelse(sf_df$Access_Metric >= cutoff, "Service desert", "Better access")
      } else {
        sf_df$Desert_Status <- ifelse(sf_df$Access_Metric <= cutoff, "Service desert", "Better access")
      }
    } else {
      if (method == "Nearest service distance") {
        cutoff <- stats::quantile(sf_df$Access_Metric, 2/3, na.rm = TRUE)
        sf_df$Desert_Status <- ifelse(sf_df$Access_Metric >= cutoff, "Service desert", "Better access")
      } else {
        cutoff <- stats::quantile(sf_df$Access_Metric, 1/3, na.rm = TRUE)
        sf_df$Desert_Status <- ifelse(sf_df$Access_Metric <= cutoff, "Service desert", "Better access")
      }
    }

    model <- tryCatch(lm(outcome ~ Access_Metric, data = sf::st_drop_geometry(sf_df)), error = function(e) NULL)
    model_tbl <- if (is.null(model)) {
      data.frame(Term = "Model failed", Estimate = NA, CI = NA, P_value = NA, Interpretation = "Not estimated")
    } else {
      ci <- suppressWarnings(confint(model))
      broom::tidy(model) %>%
        dplyr::mutate(CI = paste0("[", round(ci[,1], 3), ", ", round(ci[,2], 3), "]"),
                      P_value = sapply(p.value, format_p_apa),
                      Interpretation = sapply(p.value, interpret_p)) %>%
        dplyr::select(Term = term, Estimate = estimate, CI, P_value, Interpretation)
    }

    summary_tbl <- sf::st_drop_geometry(sf_df) %>%
      dplyr::mutate(Region = Match_Key) %>%
      dplyr::arrange(dplyr::desc(Desert_Status), Access_Metric) %>%
      dplyr::select(Region, n_obs, Outcome = outcome, Access_Metric, Desert_Status) %>%
      dplyr::mutate(Access_Metric = round(Access_Metric, 3), Outcome = round(Outcome, 4))

    desert_tbl <- sf::st_drop_geometry(sf_df) %>%
      dplyr::group_by(Desert_Status) %>%
      dplyr::summarise(Regions = dplyr::n(), Mean_Outcome = mean(outcome, na.rm = TRUE),
                       Mean_Access_Metric = mean(Access_Metric, na.rm = TRUE), .groups = "drop") %>%
      dplyr::mutate(Mean_Outcome = round(Mean_Outcome, 4), Mean_Access_Metric = round(Mean_Access_Metric, 3))

    list(data = sf_df, summary = summary_tbl, desert = desert_tbl, model = model_tbl, method = method)
  })

  output$p11_summary_table <- renderTable({ req(paper11_res()); paper11_res()$summary }, striped = TRUE, bordered = TRUE, spacing = "m")
  output$p11_desert_table <- renderTable({ req(paper11_res()); paper11_res()$desert }, striped = TRUE, bordered = TRUE)
  output$p11_model_table <- renderTable({ req(paper11_res()); paper11_res()$model }, striped = TRUE, bordered = TRUE, digits = 4)
  output$p11_map <- renderPlot({
    req(paper11_res())
    ggplot(paper11_res()$data) +
      geom_sf(aes(fill = Desert_Status), color = "white", size = 0.2) +
      scale_fill_manual(values = c("Service desert" = "#c0392b", "Better access" = "#27ae60"), na.value = "grey80") +
      theme_void() + labs(title = "Service Desert / Accessibility Classification", subtitle = paper11_res()$method, fill = "Status")
  })
  output$p11_gradient_plot <- renderPlot({
    req(paper11_res())
    df <- sf::st_drop_geometry(paper11_res()$data)
    ggplot(df, aes(x = Access_Metric, y = outcome)) +
      geom_point(aes(color = Desert_Status), size = 3, alpha = 0.8) +
      geom_smooth(method = "lm", se = TRUE, color = "black") +
      scale_color_manual(values = c("Service desert" = "#c0392b", "Better access" = "#27ae60")) +
      theme_minimal(base_size = 13) + labs(title = "Accessibility-Outcome Gradient", x = unique(df$Access_Label)[1], y = "Outcome")
  })

  paper12_res <- eventReactive(input$run_p12, {
    req(adal_spatial(), processed_df(), input$p12_predictor)
    dfp <- processed_df()
    dfp$wgt_f <- if(input$is_dhs_wgt && input$v_wgt %in% names(dfp)) dfp[[input$v_wgt]]/1000000 else if(input$v_wgt %in% names(dfp)) as.numeric(dfp[[input$v_wgt]]) else 1
    ag <- dfp %>%
      dplyr::mutate(pred_num = suppressWarnings(as.numeric(as.character(.data[[input$p12_predictor]])))) %>%
      dplyr::group_by(!!sym(input$v_reg)) %>%
      dplyr::summarise(Local_X = weighted.mean(pred_num, wgt_f, na.rm = TRUE), .groups = "drop") %>%
      dplyr::mutate(Match_Key = stringr::str_to_upper(stringr::str_trim(as.character(!!sym(input$v_reg)))))
    sf_df <- adal_spatial() %>% dplyr::left_join(ag %>% dplyr::select(Match_Key, Local_X), by = "Match_Key") %>% dplyr::filter(!is.na(outcome), !is.na(Local_X))
    validate(need(nrow(sf_df) >= 4, "At least four mapped regions with valid predictor values are required."))
    coords <- sf::st_coordinates(sf::st_centroid(sf_df))
    k <- min(input$p12_k, nrow(sf_df) - 1)
    nb <- spdep::knn2nb(spdep::knearneigh(coords, k = k))
    global_mod <- lm(outcome ~ Local_X, data = sf::st_drop_geometry(sf_df))
    local_rows <- lapply(seq_len(nrow(sf_df)), function(i) {
      idx <- unique(c(i, nb[[i]]))
      dloc <- sf::st_drop_geometry(sf_df)[idx, ]
      if (nrow(dloc) < 3 || length(unique(dloc$Local_X)) < 2) {
        return(data.frame(Region = sf_df$Match_Key[i], Local_Coefficient = NA, Local_R2 = NA, P_value = NA, Significance = "-"))
      }
      mod <- tryCatch(lm(outcome ~ Local_X, data = dloc), error = function(e) NULL)
      if (is.null(mod)) return(data.frame(Region = sf_df$Match_Key[i], Local_Coefficient = NA, Local_R2 = NA, P_value = NA, Significance = "-"))
      sm <- summary(mod)
      pval <- coef(sm)[2,4]
      data.frame(Region = sf_df$Match_Key[i], Local_Coefficient = coef(mod)[2], Local_R2 = sm$r.squared, P_value = pval, Significance = sig_label(pval))
    }) %>% dplyr::bind_rows()
    sf_df <- sf_df %>% dplyr::left_join(local_rows, by = c("Match_Key" = "Region"))
    g <- broom::tidy(global_mod, conf.int = TRUE) %>%
      dplyr::mutate(P_value = sapply(p.value, format_p_apa), Significance = sapply(p.value, sig_label)) %>%
      dplyr::select(Term = term, Estimate = estimate, CI_Lower = conf.low, CI_Upper = conf.high, P_value, Significance)
    list(data = sf_df, global = g, local = local_rows)
  })
  output$p12_global_table <- renderTable({ req(paper12_res()); paper12_res()$global }, striped = TRUE, bordered = TRUE, digits = 4)
  output$p12_local_table <- renderTable({ req(paper12_res()); paper12_res()$local %>% dplyr::mutate(P_value = sapply(P_value, format_p_apa)) }, striped = TRUE, bordered = TRUE, digits = 4)
  output$p12_coef_map <- renderPlot({ req(paper12_res()); ggplot(paper12_res()$data) + geom_sf(aes(fill = Local_Coefficient), color = "white", size = .2) + scale_fill_gradient2(low="#2980b9", mid="white", high="#c0392b") + theme_void() + labs(title="Local Coefficient Map", fill="Local beta") })
  output$p12_r2_map <- renderPlot({ req(paper12_res()); ggplot(paper12_res()$data) + geom_sf(aes(fill = Local_R2), color = "white", size = .2) + scale_fill_viridis_c(option="mako") + theme_void() + labs(title="Local Model Fit Map", fill="Local R-squared") })

  paper13_res <- eventReactive(input$run_p13, {
    req(adal_spatial())
    sf_df <- adal_spatial() %>% dplyr::filter(!is.na(outcome))
    indicator_tbl <- list()
    score_mat <- data.frame(Match_Key = sf_df$Match_Key)
    if (isTRUE(input$p13_include_outcome)) {
      score_mat$outcome_scaled <- scales::rescale(sf_df$outcome, to = c(0, 100), from = range(sf_df$outcome, na.rm = TRUE))
      indicator_tbl[[length(indicator_tbl)+1]] <- data.frame(Indicator = "Mapped outcome", Weight = 1, Mean_Scaled_Score = mean(score_mat$outcome_scaled, na.rm = TRUE))
    }
    inds <- input$p13_indicators
    if (!is.null(inds) && length(inds) > 0) {
      dfp <- processed_df()
      dfp$wgt_f <- if(input$is_dhs_wgt && input$v_wgt %in% names(dfp)) dfp[[input$v_wgt]]/1000000 else if(input$v_wgt %in% names(dfp)) as.numeric(dfp[[input$v_wgt]]) else 1
      for (v in inds) {
        ag <- dfp %>% dplyr::mutate(v_num = suppressWarnings(as.numeric(as.character(.data[[v]])))) %>%
          dplyr::group_by(!!sym(input$v_reg)) %>% dplyr::summarise(v_mean = weighted.mean(v_num, wgt_f, na.rm = TRUE), .groups = "drop") %>%
          dplyr::mutate(Match_Key = stringr::str_to_upper(stringr::str_trim(as.character(!!sym(input$v_reg)))))
        tmp <- sf_df %>% sf::st_drop_geometry() %>% dplyr::select(Match_Key) %>% dplyr::left_join(ag %>% dplyr::select(Match_Key, v_mean), by = "Match_Key")
        if (any(!is.na(tmp$v_mean)) && diff(range(tmp$v_mean, na.rm=TRUE)) > 0) {
          cname <- make.names(paste0(v, "_scaled"))
          score_mat[[cname]] <- scales::rescale(tmp$v_mean, to = c(0, 100), from = range(tmp$v_mean, na.rm = TRUE))
          indicator_tbl[[length(indicator_tbl)+1]] <- data.frame(Indicator = v, Weight = 1, Mean_Scaled_Score = mean(score_mat[[cname]], na.rm = TRUE))
        }
      }
    }
    numeric_cols <- names(score_mat)[names(score_mat) != "Match_Key"]
    validate(need(length(numeric_cols) > 0, "Select at least one indicator or include the outcome."))
    score_mat$Priority_Score <- rowMeans(score_mat[numeric_cols], na.rm = TRUE)
    score_mat$Priority_Quintile <- dplyr::ntile(score_mat$Priority_Score, 5)
    sf_df <- sf_df %>% dplyr::left_join(score_mat %>% dplyr::select(Match_Key, Priority_Score, Priority_Quintile), by = "Match_Key")
    priority_tbl <- sf::st_drop_geometry(sf_df) %>% dplyr::arrange(dplyr::desc(Priority_Score)) %>%
      dplyr::mutate(Rank = dplyr::row_number(), Priority_Score = round(Priority_Score, 2)) %>%
      dplyr::select(Rank, Region = Match_Key, n_obs, Outcome = outcome, Priority_Score, Priority_Quintile)
    indicator_df <- dplyr::bind_rows(indicator_tbl) %>% dplyr::mutate(Mean_Scaled_Score = round(Mean_Scaled_Score, 2))
    list(data = sf_df, priority = priority_tbl, indicators = indicator_df)
  })
  output$p13_priority_table <- renderTable({ req(paper13_res()); paper13_res()$priority }, striped = TRUE, bordered = TRUE, digits = 4)
  output$p13_indicator_table <- renderTable({ req(paper13_res()); paper13_res()$indicators }, striped = TRUE, bordered = TRUE)
  output$p13_priority_map <- renderPlot({ req(paper13_res()); ggplot(paper13_res()$data) + geom_sf(aes(fill = Priority_Score), color = "white", size = .2) + scale_fill_viridis_c(option="magma", direction=-1) + theme_void() + labs(title="Spatial Priority Index", fill="Score") })
  output$p13_contribution_plot <- renderPlot({ req(paper13_res()); ggplot(paper13_res()$indicators, aes(x = reorder(Indicator, Mean_Scaled_Score), y = Mean_Scaled_Score)) + geom_col(fill="#006400") + coord_flip() + theme_minimal(base_size=13) + labs(title="Indicator Contribution Summary", x="", y="Mean scaled score") })

  paper14_res <- eventReactive(input$run_p14, {
    req(adal_spatial())
    sf_df <- adal_spatial() %>% dplyr::filter(!is.na(outcome))
    nb <- spdep::poly2nb(sf_df, queen = TRUE)
    if(any(spdep::card(nb) == 0)) nb <- spdep::knn2nb(spdep::knearneigh(sf::st_coordinates(sf::st_centroid(sf_df)), k = 1))
    lw <- spdep::nb2listw(nb, style = "W", zero.policy = TRUE)
    mor <- spdep::moran.test(sf_df$outcome, lw, zero.policy = TRUE)
    gi <- spdep::localG(sf_df$outcome, lw, zero.policy = TRUE)
    loc <- spdep::localmoran(sf_df$outcome, lw, zero.policy = TRUE)
    thr <- stats::quantile(sf_df$outcome, input$p14_threshold/100, na.rm = TRUE)
    sf_df$High_Burden <- sf_df$outcome >= thr
    sf_df$GiZ <- as.numeric(gi)
    sf_df$LISA_P <- loc[,5]
    sf_df$Cluster_Class <- dplyr::case_when(
      sf_df$High_Burden & sf_df$GiZ >= 1.65 ~ "High-burden hot cluster",
      sf_df$High_Burden ~ "High-burden candidate",
      sf_df$GiZ <= -1.65 ~ "Cold cluster",
      TRUE ~ "Not clustered"
    )
    global_tbl <- data.frame(
      Test = c("Global Moran's I", "High-burden threshold", "Mean Getis-Ord Gi* z", "Number of high-burden candidate areas"),
      Statistic = c(round(as.numeric(mor$estimate[1]), 4), round(thr, 4), round(mean(sf_df$GiZ, na.rm=TRUE), 4), sum(sf_df$High_Burden, na.rm=TRUE)),
      P_value = c(format_p_apa(mor$p.value), "-", "-", "-"),
      Interpretation = c(interpret_p(mor$p.value), paste0("Top ", 100-input$p14_threshold, "% highest outcome areas"), "Average local hot/cold score", "Candidate cluster count")
    )
    local_tbl <- sf::st_drop_geometry(sf_df) %>% dplyr::arrange(dplyr::desc(outcome)) %>%
      dplyr::select(Region = Match_Key, n_obs, Outcome = outcome, GiZ, LISA_P, Cluster_Class) %>%
      dplyr::mutate(Outcome = round(Outcome, 4), GiZ = round(GiZ, 3), LISA_P = sapply(LISA_P, format_p_apa))
    profile_tbl <- sf::st_drop_geometry(sf_df) %>% dplyr::group_by(Cluster_Class) %>%
      dplyr::summarise(Regions = dplyr::n(), Mean_Outcome = mean(outcome, na.rm=TRUE), .groups="drop")
    list(data = sf_df, global = global_tbl, local = local_tbl, profile = profile_tbl)
  })
  output$p14_global_table <- renderTable({ req(paper14_res()); paper14_res()$global }, striped = TRUE, bordered = TRUE)
  output$p14_local_table <- renderTable({ req(paper14_res()); paper14_res()$local }, striped = TRUE, bordered = TRUE)
  output$p14_cluster_map <- renderPlot({ req(paper14_res()); ggplot(paper14_res()$data) + geom_sf(aes(fill = Cluster_Class), color="white", size=.2) + scale_fill_manual(values=c("High-burden hot cluster"="#b2182b", "High-burden candidate"="#ef8a62", "Cold cluster"="#2166ac", "Not clustered"="#d9d9d9"), na.value="grey80") + theme_void() + labs(title="Spatial Cluster Detection Map", fill="Cluster") })
  output$p14_cluster_plot <- renderPlot({ req(paper14_res()); ggplot(paper14_res()$profile, aes(x = reorder(Cluster_Class, Mean_Outcome), y = Mean_Outcome)) + geom_col(fill="#006400") + coord_flip() + theme_minimal(base_size=13) + labs(title="Cluster Burden Profile", x="", y="Mean outcome") })

  output$p15_vulnerable_ui <- renderUI({
    req(input$p15_equity_group, processed_df())
    choices <- unique(as.character(haven::as_factor(processed_df()[[input$p15_equity_group]])))
    choices <- choices[!is.na(choices) & choices != ""]
    selectInput("p15_vulnerable_category", "Vulnerable/reference category to map:", choices = choices)
  })

  paper15_res <- eventReactive(input$run_p15, {
    req(processed_df(), adal_spatial(), input$p15_equity_group, input$p15_vulnerable_category)
    df <- processed_df()
    df$wgt_f <- if(input$is_dhs_wgt && input$v_wgt %in% names(df)) df[[input$v_wgt]]/1000000 else if(input$v_wgt %in% names(df)) as.numeric(df[[input$v_wgt]]) else 1
    df$eq_group <- as.character(haven::as_factor(df[[input$p15_equity_group]]))
    group_tbl <- df %>% dplyr::filter(!is.na(eq_group), eq_group != "", !is.na(dv_final)) %>%
      dplyr::group_by(eq_group) %>%
      dplyr::summarise(n = dplyr::n(), Weighted_N = sum(wgt_f, na.rm=TRUE), Outcome = weighted.mean(dv_final, wgt_f, na.rm=TRUE), .groups="drop") %>%
      dplyr::arrange(dplyr::desc(Outcome))
    max_out <- max(group_tbl$Outcome, na.rm=TRUE); min_out <- min(group_tbl$Outcome, na.rm=TRUE)
    gap_tbl <- data.frame(
      Metric = c("Highest group outcome", "Lowest group outcome", "Absolute equity gap", "Relative equity ratio"),
      Estimate = c(max_out, min_out, max_out - min_out, ifelse(min_out == 0, NA, max_out/min_out))
    ) %>% dplyr::mutate(Estimate = round(Estimate, 4))
    region_tbl <- df %>% dplyr::filter(!is.na(eq_group), eq_group != "") %>%
      dplyr::mutate(vuln = eq_group == input$p15_vulnerable_category) %>%
      dplyr::group_by(!!sym(input$v_reg)) %>%
      dplyr::summarise(Vulnerable_Share = weighted.mean(as.numeric(vuln), wgt_f, na.rm=TRUE),
                       Outcome = weighted.mean(dv_final, wgt_f, na.rm=TRUE), n = dplyr::n(), .groups="drop") %>%
      dplyr::mutate(Match_Key = stringr::str_to_upper(stringr::str_trim(as.character(!!sym(input$v_reg)))))
    map_df <- adal_spatial() %>% dplyr::left_join(region_tbl %>% dplyr::select(Match_Key, Vulnerable_Share), by = "Match_Key")
    reg_out <- region_tbl %>% dplyr::arrange(dplyr::desc(Vulnerable_Share)) %>% dplyr::select(Region = Match_Key, n, Vulnerable_Share, Outcome) %>% dplyr::mutate(Vulnerable_Share = round(Vulnerable_Share,4), Outcome = round(Outcome,4))
    list(group = group_tbl, gap = gap_tbl, region = reg_out, map = map_df)
  })
  output$p15_gap_table <- renderTable({ req(paper15_res()); paper15_res()$gap }, striped = TRUE, bordered = TRUE)
  output$p15_region_table <- renderTable({ req(paper15_res()); paper15_res()$region }, striped = TRUE, bordered = TRUE)
  output$p15_share_map <- renderPlot({ req(paper15_res()); ggplot(paper15_res()$map) + geom_sf(aes(fill = Vulnerable_Share), color="white", size=.2) + scale_fill_viridis_c(option="magma", labels=scales::percent) + theme_void() + labs(title="Vulnerable Group Share by Region", fill="Share") })
  output$p15_group_plot <- renderPlot({ req(paper15_res()); ggplot(paper15_res()$group, aes(x = reorder(eq_group, Outcome), y = Outcome)) + geom_col(fill="#006400") + coord_flip() + theme_minimal(base_size=13) + labs(title="Outcome by Equity Group", x="", y="Weighted outcome") })


  #########################################################################################################################
  ###############################################################################
  # --- FINALIZED REPORTING ENGINE (Fixed YAML Params) ---
  output$dl_report <- downloadHandler(
    filename = function() {
      paste0("HawaSpatial_Full_Report_", Sys.Date(), ".", input$format)
    },
    content = function(file) {
      # Use a notification so the user knows it's working
      id <- showNotification("Generating Report... This may take a moment.", duration = NULL, closeButton = FALSE)
      on.exit(removeNotification(id), add = TRUE)

      tempReport <- file.path(tempdir(), "hawa_final_report.Rmd")

      # Helper to safely grab reactive values for the report
      safe_get <- function(expr) {
        tryCatch({ expr }, error = function(e) { NULL })
      }
      # We construct the RMarkdown file line-by-line as a character vector
      report_content <- c(
        "---",
        "title: '`r params$doc_title`'",
        "author: 'HawaSpatial Pro Intelligence Platform'",
        "date: '`r Sys.Date()`'",
        "output:",
        "  word_document: ",
        "    fig_caption: yes",
        "    fig_width: 9",
        "    fig_height: 6",
        "  pdf_document: default",
        "  html_document: default",
        "params:",
        "  doc_title: NA",
        "  spatial_data: NA",
        "  models: NA",
        "  outcome_type: NA",
        "  v_reg: NA",
        "  v_cluster: NA",
        "  clean_df: NA",
        "  st_data: NA",
        "  v_time: NA",
        "  v_wealth: NA",
        "  raw_df: NA",
        "  v_indiv_cont: NA",
        "  v_indiv_fact: NA",
        "  v_comm_cont: NA",
        "  v_comm_fact: NA",
        "  v_wgt: NA",
        "  is_dhs_wgt: NA",
        "  cov_vars: NA",
        "  sae_data: NA",
        "  econ_data: NA",
        "  econ_type: NA",
        "  econ_lm: NA",
        "  econ_fit_table: NA",
        "  econ_weights_table: NA",
        "  econ_resid_table: NA",
        "  gini_val: NA",
        "  theil_val: NA",
        "  outcome_vec: NA",
        "  decomp_data: NA",
        "  spatial_decomp_data: NA",
        "  conc_val: NA",
        "  joint_model: NA",
        "  joint_spatial: NA",
        "  joint_decomp_data: NA",
        "  joint_st_decomp_data: NA",
        "  st_int_model: NA",
        "  st_int_data: NA",
        "  st_moran_table: NA",
        "  st_rate_table: NA",
        "  st_lisa_table: NA",
        "  st_hotspot_table: NA",
        "  bivariate_table: NA",
        "  ineq_apa_table: NA",
        "  ineq_region_table: NA",
        "  paper11_data: NA",
        "  paper12_data: NA",
        "  paper13_data: NA",
        "  paper14_data: NA",
        "  paper15_data: NA",
        # ---------------------------------
        "---",
        "",
        "```{r setup, include=FALSE}",
        "knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)",
        "# Load necessary libraries inside the report environment",
        "library(tidyverse); library(sf); library(spdep); library(glmmTMB); library(broom.mixed)",
        "library(knitr); library(cowplot); library(biscale); library(gstat); library(spatialreg); library(classInt); library(scales); library(corrplot); library(RColorBrewer)",
        "",
        "# Define ICC function for report usage (MATCHING SHINY UI LOGIC)",
        "calc_icc_internal <- function(model) {",
        "  tryCatch({",
        "    vc <- VarCorr(model)",
        "    var_u <- as.numeric(vc$cond$cluster_idx[1,1])",
        "    fam <- family(model)$family",
        "    if (grepl('binomial', fam)) { var_e <- 3.289868 }",
        "    else if (grepl('poisson', fam) || grepl('nbinom', fam)) { var_e <- 1.0 }",
        "    else if (grepl('beta', fam)) { var_e <- 1.0 }",
        "    else { var_e <- sigma(model)^2 }",
        "    icc_val <- var_u / (var_u + var_e)",
        "    return(icc_val)",
        "  }, error = function(e) return(0))",
        "}",
        "# Define PCV function for report usage",
        "calc_pcv_internal <- function(m0, m_full) {",
        "  icc_m0 <- calc_icc_internal(m0)",
        "  icc_m_full <- calc_icc_internal(m_full)",
        "  if (icc_m0 > 0) {",
        "    pcv_val <- (icc_m0 - icc_m_full) / icc_m0",
        "    return(max(0, pcv_val))",
        "  } else {",
        "    return(0)",
        "  }",
        "}",
        "shp_geom <- params$spatial_data %>% dplyr::select(Match_Key, geometry) %>% distinct()",
        "```",
        "",
        "# 1. Introduction and 15-Paper Publication Roadmap",
        "HawaSpatial Pro is structured as an integrated spatial-analysis publication engine. After this introduction, the report presents the complete paper suite in sequential order, from Paper 1 to Paper 15, while also retaining the foundational descriptive, ESDA, bivariate, advanced geostatistical, and model-diagnostic outputs needed to support each manuscript.",
        "",
        "The report covers: Paper 1 Small Area Estimation; Paper 2 Spatial Econometrics; Paper 3 Spatial Inequality; Paper 4 Multilevel Determinants; Paper 5 Regional Decomposition; Paper 6 Spatiotemporal Dynamics and Emerging Hotspots; Paper 7 Spatiotemporal Decomposition; Paper 8 Joint Spatial Modeling; Paper 9 Joint Regional Decomposition; Paper 10 Joint Spatiotemporal Decomposition; Paper 11 Spatial Accessibility and Service Deserts; Paper 12 Local Determinants; Paper 13 Spatial Priority Index; Paper 14 Spatial Scan and Cluster Detection; and Paper 15 Spatial Equity and Environmental Justice.",
        "",
        "## 1.1 Foundational Descriptive Table: Characteristics of Respondents",
        "Descriptive statistics (Weighted Frequency/Percentage for categorical, Weighted Mean/SD for continuous).",
        "",
        "```{r table1_desc, results='asis'}",
        "# Use RAW data to get the labels (Urban/Rural) instead of numbers",
        "df_rep <- params$raw_df",
        "wt_rep <- if(params$is_dhs_wgt) df_rep[[params$v_wgt]]/1000000 else df_rep[[params$v_wgt]]",
        "if(is.null(wt_rep)) wt_rep <- rep(1, nrow(df_rep))",
        "table_rows <- list()",
        "",
        "# Continuous",
        "cont_vars <- c(params$v_indiv_cont, params$v_comm_cont)",
        "for(v in cont_vars) {",
        "  val <- as.numeric(df_rep[[v]]); valid <- !is.na(val) & !is.na(wt_rep)",
        "  if(sum(valid) > 0) {",
        "    w_mean <- weighted.mean(val[valid], wt_rep[valid])",
        "    w_sd <- sqrt(cov.wt(data.frame(val[valid]), wt = wt_rep[valid])$cov[1,1])",
        "    table_rows[[length(table_rows)+1]] <- data.frame(Variable=v, Category='Mean (SD)', Statistic=paste0(round(w_mean,2), ' (', round(w_sd,2), ')'))",
        "  }",
        "}",
        "",
        "# Categorical",
        "cat_vars <- c(params$v_indiv_fact, params$v_comm_fact)",
        "for(v in cat_vars) {",
        "  # Convert to factor to get labels",
        "  val_fac <- haven::as_factor(df_rep[[v]])",
        "  tmp <- data.frame(cat = val_fac, w = wt_rep) %>% drop_na()",
        "  if(nrow(tmp) > 0) {",
        "     agg <- tmp %>% group_by(cat) %>% summarise(n_w = sum(w)) %>% mutate(pct = n_w/sum(n_w)*100)",
        "     table_rows[[length(table_rows)+1]] <- data.frame(Variable=paste0('**', v, '**'), Category='', Statistic='')",
        "     for(i in 1:nrow(agg)) {",
        "       table_rows[[length(table_rows)+1]] <- data.frame(Variable='', Category=as.character(agg$cat[i]), Statistic=paste0(round(agg$n_w[i],0), ' (', round(agg$pct[i],1), '%)'))",
        "     }",
        "  }",
        "}",
        "",
        "if(length(table_rows) > 0) {",
        "  final_t1 <- do.call(rbind, table_rows)",
        "  kable(final_t1, caption = 'Weighted Descriptive Statistics', align = 'l')",
        "} else { cat('No variables selected for Table 1.') }",
        "```",
        "## 3. Methodological Transparency",
        "All analyses were performed using the HawaSpatial Pro Intelligence Platform. The underlying statistical engines utilize the following R packages:",
        "* **Multilevel Modeling:** `glmmTMB` for hierarchical estimation.",
        "* **Spatial Econometrics:** `spatialreg` and `spdep` for SAR/SEM and Lagrange Multiplier tests.",
        "* **Small Area Estimation:** `sae` for Fay-Herriot area-level models.",
        "* **Spatial Inequality:** `ineq` for Gini and Theil index calculations.",
        "* **Survey Weighting:** `survey` for design-based inference.",
        "## 4. Exploratory Spatial Data Analysis (ESDA)",
        "This section reproduces all 17 spatial diagnostic methods.",
        "",
        "```{r esda_loop, results='asis'}",
        "df <- params$spatial_data",
        "nb <- poly2nb(df, queen = TRUE)",
        "if(any(card(nb) == 0)) nb <- knn2nb(knearneigh(st_coordinates(st_centroid(df)), k = 1))",
        "lw <- nb2listw(nb, style = 'W', zero.policy = TRUE)",
        "df_norm <- data.frame(x = seq(-4, 4, length.out = 1000)) %>% mutate(y = dnorm(x))",
        "",
        "methods <- c('1. Connectivity Graph', '2. Basic Choropleth', '3. EBS Smoothing', ",
        "             '4. Quantile Map', '5. Global Moran\\'s I', '6. Moran Scatterplot', ",
        "             '7. LISA Cluster Map', '8. Getis-Ord Gi*', '9. Join Count Analysis', ",
        "             '10. Hinge Map', '11. Spatial Correlogram', '12. Regime Density', ",
        "             '13. Spatial Lag Map', '14. SDG Target Gap', '15. Local Geary\\'s C', ",
        "             '16. Conditional Plot', '17. Global Geary\\'s C')",
        "",
        "for (m in methods) {",
        "  cat('\\n\\n## ', m, '\\n')",
        "  ",
        "  tryCatch({",
        "    if(m == '1. Connectivity Graph') {",
        "      plot(st_geometry(df), border='grey')",
        "      plot(nb, st_coordinates(st_centroid(df)), add=TRUE, col='#006400', lwd=2)",
        "    } else if(m == '2. Basic Choropleth') {",
        "      p <- ggplot(df) + geom_sf(aes(fill=outcome)) + scale_fill_distiller(palette = 'RdPu', direction = 1, labels = percent) + theme_void()",
        "      print(p)",
        "    } else if(m == '3. EBS Smoothing') {",
        "      eb_res <- EBlocal(ri = df$cases, ni = df$n_obs, nb = nb, zero.policy = TRUE)",
        "      df$eb_rate <- eb_res$est",
        "      p1 <- ggplot(df) + geom_sf(aes(fill=outcome)) + ",
        "        geom_sf_text(aes(label = Match_Key), size = 2, check_overlap = TRUE) + ",
        "        scale_fill_distiller(palette='RdPu', direction=1, labels = percent) + ",
        "        labs(title='Raw Rates') + theme_void() + theme(legend.position='bottom')",
        "      p2 <- ggplot(df) + geom_sf(aes(fill=eb_rate)) + ",
        "        geom_sf_text(aes(label = Match_Key), size = 2, check_overlap = TRUE) + ",
        "        scale_fill_distiller(palette='RdPu', direction=1, labels = percent) + ",
        "        labs(title='EBS Smoothed') + theme_void() + theme(legend.position='bottom')",
        "      print(plot_grid(p1, p2, ncol = 2))",
        "    } else if(m == '4. Quantile Map') {",
        "      breaks_qt <- classIntervals(df$outcome, n = 5, style = 'quantile')",
        "      df$quantile_cat <- cut(df$outcome, breaks_qt$brks, include.lowest = TRUE)",
        "      p <- ggplot(df) + geom_sf(aes(fill = quantile_cat)) + scale_fill_brewer(palette = 'Blues') + theme_void()",
        "      print(p)",
        "    } else if(m == '5. Global Moran\\'s I') {",
        "      res <- moran.test(df$outcome, lw, zero.policy = TRUE)",
        "      z <- as.numeric(res$statistic)",
        "      p <- ggplot(df_norm, aes(x, y)) + geom_line(size = 1) + ",
        "        geom_area(data = subset(df_norm, x >= z), fill = 'red', alpha = 0.5) +",
        "        geom_vline(xintercept = z, color = 'red', linetype='dashed') +",
        "        annotate('text', x = 2.5, y = 0.3, label = paste('p:', format.pval(res$p.value)), color='red', fontface='bold') +",
        "        labs(subtitle = paste('Moran\\'s I:', round(res$estimate[1],3), '| Z:', round(z,2))) + theme_minimal()",
        "      print(p)",
        "    } else if(m == '6. Moran Scatterplot') {",
        "      df$lag_outcome <- lag.listw(lw, df$outcome, zero.policy = TRUE)",
        "      m_prev <- mean(df$outcome)",
        "      df$quadrant <- case_when(df$outcome >= m_prev & df$lag_outcome >= m_prev ~ 'High-High', df$outcome <= m_prev & df$lag_outcome <= m_prev ~ 'Low-Low', df$outcome >= m_prev & df$lag_outcome <= m_prev ~ 'High-Low', df$outcome <= m_prev & df$lag_outcome >= m_prev ~ 'Low-High')",
        "      p <- ggplot(df, aes(x = outcome, y = lag_outcome, color = quadrant)) + ",
        "        geom_vline(xintercept = m_prev, linetype = 'dashed') + geom_hline(yintercept = mean(df$lag_outcome), linetype = 'dashed') + ",
        "        geom_point(size = 3) + geom_text(aes(label = Match_Key), vjust = -0.5, size = 2.5, check_overlap = TRUE, show.legend = FALSE) + ",
        "        scale_color_manual(values = c('High-High'='red', 'Low-Low'='blue', 'High-Low'='pink', 'Low-High'='lightblue')) + theme_minimal()",
        "      print(p)",
        "    } else if(m == '7. LISA Cluster Map') {",
        "      loc_m <- localmoran(df$outcome, lw, zero.policy = TRUE)",
        "      m_prev <- mean(df$outcome)",
        "      df$lag_outcome <- lag.listw(lw, df$outcome, zero.policy = TRUE)",
        "      df$Pr <- loc_m[,5]",
        "      df$lisa_quadrant <- case_when(",
        "        df$Pr <= 0.05 & df$outcome >= m_prev & df$lag_outcome >= m_prev ~ 'High-High', ",
        "        df$Pr <= 0.05 & df$outcome <= m_prev & df$lag_outcome <= m_prev ~ 'Low-Low', ",
        "        df$Pr <= 0.05 & df$outcome >= m_prev & df$lag_outcome <= m_prev ~ 'High-Low',",
        "        df$Pr <= 0.05 & df$outcome <= m_prev & df$lag_outcome >= m_prev ~ 'Low-High',",
        "        TRUE ~ 'Not Significant'",
        "      )",
        "      lisa_colors <- c('High-High'='red', 'Low-Low'='blue', 'High-Low'='pink', 'Low-High'='lightblue', 'Not Significant'='lightgray')",
        "      p <- ggplot(df) + geom_sf(aes(fill=lisa_quadrant), color='black') + ",
        "        scale_fill_manual(values=lisa_colors, name='LISA Cluster Type') + ",
        "        labs(title = '7. LISA Cluster Map (Local Spatial Association)') + theme_void()",
        "      print(p)",
        "    } else if(m == '8. Getis-Ord Gi*') {",
        "      gi_stat <- localG(df$outcome, lw, zero.policy = TRUE)",
        "      p <- ggplot(df) + geom_sf(aes(fill=as.numeric(gi_stat))) + scale_fill_gradient2() + theme_void()",
        "      print(p)",
        "    } else if(m == '9. Join Count Analysis') {",
        "      df$Coverage_Class <- factor(ifelse(df$outcome > mean(df$outcome), 'High', 'Low'))",
        "      p <- ggplot(df) + geom_sf(aes(fill = Coverage_Class)) + scale_fill_manual(values = c('High' = 'forestgreen', 'Low' = 'darkorange')) + theme_void()",
        "      print(p)",
        "    } else if(m == '10. Hinge Map') {",
        "      cuts <- quantile(df$outcome, c(0, 0.25, 0.5, 0.75, 1)); iqr <- diff(cuts[c(2, 4)])",
        "      df$Hinge <- cut(df$outcome, breaks = c(-Inf, cuts[2]-1.5*iqr, cuts[2], cuts[3], cuts[4], cuts[4]+1.5*iqr, Inf), labels = c('L-Outlier', '<25%', '25-50%', '50-75%', '>75%', 'U-Outlier'))",
        "      p <- ggplot(df) + geom_sf(aes(fill=Hinge)) + scale_fill_brewer(palette='RdBu') + theme_void()",
        "      print(p)",
        "    } else if(m == '11. Spatial Correlogram') {",
        "      n_regions <- nrow(df)",
        "      max_lag <- ifelse(n_regions < 8, floor(n_regions / 2), 4)",
        "      if(n_regions < 4) { cat('Not enough regions.') } else {",
        "         cor_res <- sp.correlogram(nb, df$outcome, order = max_lag, method = 'I', style = 'W', zero.policy = TRUE)",
        "         actual_lags <- nrow(cor_res$res)",
        "         cor_df <- data.frame(",
        "           Lag = 1:actual_lags,",
        "           Morans_I = cor_res$res[, 1],",
        "           P_Value = 2 * (1 - pnorm(abs(cor_res$res[, 1] - cor_res$res[, 2]) / cor_res$res[, 3]))",
        "         ) %>% mutate(Significance = ifelse(P_Value < 0.05, 'Significant', 'Not Significant'))",
        "         p <- ggplot(cor_df, aes(x = Lag, y = Morans_I)) +",
        "           geom_hline(yintercept = 0, linetype = 'dashed', color = 'gray') +",
        "           geom_line() +",
        "           geom_point(aes(color = Significance), size = 4) +",
        "           scale_color_manual(values = c('Significant' = 'red', 'Not Significant' = 'blue')) +",
        "           scale_x_continuous(breaks = 1:actual_lags) +",
        "           labs(title = '11. Spatial Correlogram', y = 'Moran\\'s I') +",
        "           theme_minimal()",
        "         print(p)",
        "      }",
        "    } else if(m == '12. Regime Density') {",
        "      gi_stat <- localG(df$outcome, lw, zero.policy = TRUE)",
        "      df$Regime <- ifelse(as.numeric(gi_stat) > 1.65, 'Hot Spot', ifelse(as.numeric(gi_stat) < -1.65, 'Cold Spot', 'Neutral'))",
        "      p <- ggplot(df, aes(x=outcome, fill=Regime)) + geom_density(alpha=0.5) + theme_minimal()",
        "      print(p)",
        "    } else if(m == '13. Spatial Lag Map') {",
        "      df$lag_outcome <- lag.listw(lw, df$outcome, zero.policy = TRUE)",
        "      p <- ggplot(df) + geom_sf(aes(fill = lag_outcome)) + scale_fill_distiller(palette = 'YlGnBu', direction = 1, labels = percent) + labs(title = '13. Spatial Lag Map (Neighborhood Context)', subtitle = 'Average Rate of Neighbors') + theme_void()",
        "      print(p)",
        "    } else if(m == '14. SDG Target Gap') {",
        "      df$gap <- df$outcome - 0",
        "      p <- ggplot(df) + geom_sf(aes(fill=gap)) + scale_fill_viridis_c(option='magma', direction=-1, name='Gap') + labs(subtitle = 'Distance to Zero Target') + theme_void()",
        "      print(p)",
        "    } else if(m == '15. Local Geary\\'s C') {",
        "      df$geary_c <- as.numeric(localC(df$outcome, lw))",
        "      p <- ggplot(df) + geom_sf(aes(fill = geary_c)) + scale_fill_viridis_c(option = 'magma', direction = -1) + theme_void()",
        "      print(p)",
        "    } else if(m == '16. Conditional Plot') {",
        "      df$lag_bin <- cut(lag.listw(lw, df$outcome), breaks=3, labels=c('Low Nb', 'Med Nb', 'High Nb'))",
        "      p <- ggplot(df, aes(x=lag_bin, y=outcome, fill=lag_bin)) + geom_boxplot() + theme_minimal()",
        "      print(p)",
        "    } else if(m == '17. Global Geary\\'s C') {",
        "      res <- geary.test(df$outcome, lw, zero.policy = TRUE)",
        "      z <- as.numeric(res$statistic)",
        "      p <- ggplot(df_norm, aes(x, y)) + geom_line() + ",
        "        geom_area(data = subset(df_norm, x <= z), fill = '#006400', alpha = 0.5) +",
        "        geom_vline(xintercept = z, color = 'red', linetype='dashed') +",
        "        annotate('text', x = -2.5, y = 0.3, label = paste('p:', format.pval(res$p.value)), color='#006400', fontface='bold') +",
        "        labs(subtitle = paste('Geary\\'s C:', round(res$estimate[1],3), '| Z:', round(z,2))) + theme_minimal()",
        "      print(p)",
        "    }",
        "  }, error = function(e) {",
        "    cat(paste('Plot skipped for', m, ': Insufficient data.'))",
        "  })",
        "}",
        "```",
        "",
        "\\newpage",
        "# 5. Bivariate Spatial Analysis",
        "## 5.1 Bivariate Global Moran's I Test",
        "```{r bi_moran_test}",
        "if (nrow(df) > 1 && length(lw$neighbours) > 1 && sum(card(nb)) > 0) {",
        "  # Standardize variables for Bivariate I calculation",
        "  z_y <- scale(df$outcome)[, 1]",
        "  z_x <- scale(df$bi_var)[, 1]",
        "  # Run Moran Monte Carlo test (99 simulations)",
        "  bi_moran_mc <- moran.mc(z_y * z_x, lw, nsim = 99, zero.policy = TRUE)",
        "  ",
        "  I_B_observed <- bi_moran_mc$statistic",
        "  p_value <- bi_moran_mc$p.value",
        "  ",
        "  cat('Bivariate Moran\\'s I (Outcome vs. Lagged ', params$v_wealth, ') - Permutation Test (999 simulations):\\n')",
        "  cat(paste('Observed Bivariate I:', round(I_B_observed, 4), '\\n'))",
        "  cat(paste('P-value (Simulated):', format.pval(p_value), '\\n'))",
        "  cat(paste('Alternative hypothesis: greater\\n'))",
        "  cat('\\nInterpretation: Measures if the outcome in a region is correlated with the X variable (e.g., Wealth) in neighboring regions. A significant positive I indicates spatial covariance between the two variables.')",
        "} else {",
        "  cat('Skipped: Requires at least two connected regions for Bivariate Moran\\'s I.')",
        "}",
        "```",
        "",
        "## 5.2 Bivariate Map (Outcome vs. X)",
        "```{r bivariate_map}",
        "bi_data <- bi_class(df, x = bi_var, y = outcome, style = 'quantile', dim = 3)",
        "map <- ggplot(bi_data) + geom_sf(aes(fill = bi_class), show.legend = FALSE) +",
        "  bi_scale_fill(pal = 'GrPink', dim = 3) + theme_void()",
        "legend <- bi_legend(pal = 'GrPink', dim = 3, xlab = 'Higher X', ylab = 'Higher Outcome')",
        "ggdraw() + draw_plot(map, 0, 0, 1, 1) + draw_plot(legend, 0.05, 0.05, 0.25, 0.25)",
        "```",
        "",
        "## 5.3 Covariate Correlation Matrix",
        "This plot visualizes the correlation between all independent variables selected for the Multilevel Model to check for multicollinearity.",
        "```{r cov_corr_rep, fig.height=8, fig.width=8}",
        "vars <- params$cov_vars",
        "if(length(vars) >= 2) {",
        "  # Prepare data from the clean dataframe passed in params",
        "  df_corr <- params$clean_df %>%",
        " dplyr::select(all_of(vars)) %>%",
        "    mutate(across(everything(), ~as.numeric(as.factor(.)))) %>%",
        "    na.omit()",
        "  ",
        "  M <- cor(df_corr)",
        "  ",
        "  corrplot::corrplot(M, method = 'circle', type = 'upper', ",
        "           order = 'hclust', tl.col = 'black', tl.srt = 45,",
        "           addCoef.col = 'black', number.cex = 0.6,",
        "           col = RColorBrewer::brewer.pal(n = 8, name = 'RdYlBu'))",
        "} else {",
        "  cat('Insufficient covariates selected for correlation analysis.')",
        "}",
        "```",
        "",
        "\\newpage",
        "# 6. Advanced Geostatistics",
        "## 6.1 Ordinary Kriging (Surface Interpolation)",
        "```{r krig_plot}",
        "pts <- st_centroid(df)",
        "# Handle duplicate locations for Kriging",
        "pts <- pts[!duplicated(st_coordinates(pts)), ]",
        "grid <- st_make_grid(df, n = c(50, 50), what = 'centers') %>% st_as_sf() %>% st_filter(df)",
        "v_fit <- variogram(outcome ~ 1, pts)",
        "v_mod <- fit.variogram(v_fit, model = vgm(psill=var(pts$outcome, na.rm=T), 'Exp', range=500000, nugget=0.01))",
        "krig_res <- krige(outcome ~ 1, pts, grid, model = v_mod)",
        "p_krig <- ggplot(krig_res) + geom_sf(aes(color=var1.pred), size=1) + scale_color_viridis_c() + labs(title='Kriging Surface') + theme_void()",
        "print(p_krig)",
        "```",
        "",
        "## 6.2 Spatial Lag Model Residuals Map",
        "```{r resid_plot}",
        "if (nrow(df) > 1 && length(lw$neighbours) > 1 && sum(card(nb)) > 0) {",
        "  m_sar <- lagsarlm(outcome ~ 1, data = df, lw, zero.policy = TRUE)",
        "  df$resids <- residuals(m_sar)",
        "  p_resid <- ggplot(df) + geom_sf(aes(fill=resids)) + scale_fill_gradient2() + labs(title='Spatial Lag Model Residuals') + theme_void()",
        "  print(p_resid)",
        "} else {",
        "  cat('Skipped: Requires connected regions to run Spatial Lag Model.')",
        "}",
        "```",
        "",
        "## 6.3 APA Table: Bivariate Screening Results",
        "The table below provides automatic bivariate tests before advanced spatial modeling. Categorical or factor predictors are screened using Pearson chi-square tests. Numerical predictors are screened using one-way ANOVA across outcome groups.",
        "```{r advanced_bivariate_table_rep, echo=FALSE}",
        "if(!is.null(params$bivariate_table)) {",
        "  knitr::kable(params$bivariate_table, caption = 'APA Table: Bivariate Screening Results for Advanced Spatial Analysis')",
        "} else { cat('Bivariate screening table was not available.') }",
        "```",
        "",
        "\\newpage",
        "# PAPER 1: Small Area Estimation (SAE) and Multivariate Spatial Econometrics",
        "# 7. Paper: Small Area Estimation (SAE) Analysis",
        "## 7.1 Methodology",
        "To address the limitations of direct survey estimates in regions with small sample sizes, we implemented a Fay-Herriot (FH) area-level model. This approach stabilizes regional estimates by borrowing strength from a regional covariate and accounting for the design-based sampling variance. The resulting Empirical Best Linear Unbiased Predictors (EBLUPs) provide a more reliable basis for sub-national monitoring.",
        "",
        "## 7.2 Model Diagnostics",
        "The following table presents the goodness-of-fit metrics for the Fay-Herriot model. Lower AIC and BIC values indicate a better fit between the regional covariates and the survey outcomes.",
        "```{r sae_diags_rep, echo=FALSE}",
        "if(!is.null(params$sae_data)) {",
        "  knitr::kable(params$sae_data$diags, caption = 'Fay-Herriot Model Diagnostics')",
        "} else { cat('SAE model diagnostics not available.') }",
        "```",
        "",
        "## 7.3 Spatial Distribution of Stabilized Estimates (EBLUP)",
        "```{r sae_report_plot, echo=FALSE}",
        "if(!is.null(params$sae_data)) {",
        "  df_sae_raw <- params$sae_data$data",
        "  df_sae_map <- shp_geom %>% left_join(df_sae_raw, by = 'Match_Key')",
        "  ggplot(df_sae_map) + ",
        "    geom_sf(aes(fill=eblup), color = 'white', size = 0.1) + ",
        "    scale_fill_viridis_c(option='mako', labels = scales::percent) + ",
        "    theme_void() + ",
        "    labs(title='SAE Stabilized Map (EBLUP)', subtitle='Stabilized regional prevalence rates')",
        "} else { cat('SAE model was not run.') }",
        "```",
        "",
        "## 7.4 Shrinkage Analysis: Direct vs. Model-Based Estimates",
        "The shrinkage plot visualizes how the Fay-Herriot model 'pulls' noisy direct survey estimates toward the regression line. Regions with high sampling variance experience more significant stabilization (shrinkage).",
        "```{r sae_shrink_rep, echo=FALSE, fig.height=8}",
        "if(!is.null(params$sae_data)) {",
        "  df_s <- params$sae_data$data",
        "  ggplot(df_s) +",
        "    geom_segment(aes(x = outcome, xend = eblup, y = reorder(Match_Key, outcome), yend = Match_Key), arrow = arrow(length = unit(0.15, 'cm')), color = 'gray') +",
        "    geom_point(aes(x = outcome, y = Match_Key, color = 'Direct Survey'), size = 2) +",
        "    geom_point(aes(x = eblup, y = Match_Key, color = 'SAE EBLUP'), size = 2) +",
        "    scale_color_manual(values = c('Direct Survey' = '#e74c3c', 'SAE EBLUP' = '#2980b9'), name = 'Estimate Type') +",
        "    labs(title = 'Shrinkage Effect', subtitle = 'Arrows indicate shift from direct survey data to stabilized model estimates', x = 'Prevalence', y = '') + ",
        "    theme_minimal(base_size = 10)",
        "} else { cat('SAE data not available for shrinkage plot.') }",
        "```",
        "",
        "## 7.5 Comparison of Direct and SAE Estimates",
        "The table below compares the direct survey estimates with the SAE EBLUPs, including the Coefficient of Variation (CV). A reduction in CV indicates improved precision.",
        "```{r sae_report_table, echo=FALSE}",
        "if(!is.null(params$sae_data)) {",
        "  knitr::kable(params$sae_data$data %>% dplyr::select(Region = Match_Key, `Direct Mean` = outcome, `SAE (EBLUP)` = eblup, `Direct CV (%)` = cv_direct, `SAE CV (%)` = cv_eblup), ",
        "               digits = 4, caption = 'Comparison of Precision: Direct vs. SAE')",
        "}",
        "```",
        "\\newpage",
        "# 8. Paper 2: Multivariate Spatial Econometrics & Spillover Analysis",
        "## 8.1 Methodology",
        "To account for spatial dependence, we employed Multivariate Spatial Econometric models. This approach identifies whether development gaps are localized or systemic across administrative boundaries.",
        "",
        "## 8.2 Model Selection Diagnostics (LM Tests)",
        "```{r econ_lm_report, echo=FALSE}",
        "if(!is.null(params$econ_lm)) {",
        "  lm_df <- data.frame(",
        "    Test = names(params$econ_lm),",
        "    Statistic = sapply(params$econ_lm, function(x) round(x$statistic, 3)),",
        "    P_Value = sapply(params$econ_lm, function(x) round(x$p.value, 4))",
        "  )",
        "  knitr::kable(lm_df, caption = 'Lagrange Multiplier Diagnostics')",
        "} else { cat('LM tests were not performed.') }",
        "```",
        "",
        "## 8.3 Multivariate Regression Coefficients",
        "```{r econ_coef_report, echo=FALSE}",
        "if(!is.null(params$econ_data)) {",
        "  coef_df <- as.data.frame(broom::tidy(params$econ_data$model))",
        "  knitr::kable(coef_df, digits = 4, caption = paste('Regression Results:', params$econ_type))",
        "} else { cat('Spatial regression model was not run.') }",
        "```",
        "",
        "## 8.4 Partitioning of Spillover Effects (Impacts)",
        "```{r econ_imp_report, echo=FALSE}",
        "if(!is.null(params$econ_data$impacts)) {",
        "  imp <- params$econ_data$impacts",
        "  # Extract numeric values safely",
        "  dir_v <- as.numeric(imp$direct)",
        "  ind_v <- as.numeric(imp$indirect)",
        "  tot_v <- as.numeric(imp$total)",
        "  ",
        "  # Attempt to get variable names from the impacts object or the model params",
        "  v_names <- names(imp$direct)",
        "  if(is.null(v_names)) v_names <- params$econ_data$pred_names",
        "  ",
        "  # Verify lengths match before creating data.frame to prevent the 0, 8 error",
        "  if(length(v_names) == length(dir_v) && length(dir_v) > 0) {",
        "    imp_df <- data.frame(",
        "      Variable = v_names,",
        "      Direct = dir_v,",
        "      Indirect = ind_v,",
        "      Total = tot_v",
        "    )",
        "    knitr::kable(imp_df, digits = 4, caption = 'Direct, Indirect (Spillover), and Total Effects')",
        "  } else { ",
        "    cat('Impact partitioning values were found but variable names could not be aligned.') ",
        "  }",
        "} else { cat('Impact partitioning is only available for Spatial Lag (SAR) models.') }",
        "```",
        "",
        "## 8.5 APA Model Fit and Spatial Diagnostic Table",
        "```{r econ_extra_fit_report, echo=FALSE}",
        "if(!is.null(params$econ_fit_table)) { knitr::kable(params$econ_fit_table, caption = 'APA Table: Spatial Econometric Model Fit and Diagnostics') } else { cat('Spatial econometric diagnostic table not available.') }",
        "```",
        "",
        "## 8.6 Spatial Weights Summary",
        "```{r econ_weights_report, echo=FALSE}",
        "if(!is.null(params$econ_weights_table)) { knitr::kable(params$econ_weights_table, caption = 'APA Table: Spatial Weights Matrix Summary') } else { cat('Spatial weights summary not available.') }",
        "```",
        "",
        "## 8.7 Regional Observed, Fitted, and Residual Table",
        "```{r econ_resid_report, echo=FALSE}",
        "if(!is.null(params$econ_resid_table)) { knitr::kable(params$econ_resid_table, digits = 4, caption = 'APA Table: Observed, Fitted, and Residual Values by Region') } else { cat('Spatial econometric residual table not available.') }",
        "```",
        "",
        "\\newpage",
        "# 9. Paper 3: Spatial Inequality & SDG 10 Monitoring",
        "## 9.1 Geographic Concentration Metrics",
        "Sustainable Development Goal (SDG) 10 aims to reduce inequality. We measure this via the Gini, Theil, and Concentration indices.",
        "```{r inequality_report, echo=FALSE}",
        "if(!is.null(params$gini_val)) {",
        "  ineq_df <- data.frame(",
        "    Metric = c('Gini Coefficient', 'Theil Index', 'Concentration Index'),",
        "    Value = c(round(params$gini_val, 4), round(params$theil_val, 4), round(params$conc_val, 4))",
        "  )",
        "  knitr::kable(ineq_df, caption = 'Summary of Spatial Inequality Indices')",
        "} else { cat('Inequality metrics were not calculated.') }",
        "```",
        "",
        "## 9.2 APA Table: Detailed Spatial Inequality Diagnostics",
        "```{r inequality_apa_detail_report, echo=FALSE}",
        "if(!is.null(params$ineq_apa_table)) {",
        "  params$ineq_apa_table %>% dplyr::select(Metric, Estimate = APA_Formatted, Interpretation) %>%",
        "    knitr::kable(caption = 'APA Table: Detailed Spatial Inequality Diagnostics')",
        "} else { cat('Detailed inequality diagnostics were not available.') }",
        "```",
        "",
        "## 9.3 Regional Inequality Ranking Table",
        "```{r inequality_region_table_report, echo=FALSE}",
        "if(!is.null(params$ineq_region_table)) {",
        "  params$ineq_region_table %>%",
        "    dplyr::mutate(Outcome = round(Outcome, 4), Outcome_Share_Percent = round(Outcome_Share_Percent, 2), Cumulative_Share_Percent = round(Cumulative_Share_Percent, 2), Gap_From_Mean = round(Gap_From_Mean, 4), Relative_Gap_Percent = round(Relative_Gap_Percent, 2)) %>%",
        "    knitr::kable(caption = 'APA Table: Regional Ranking and Inequality Contribution')",
        "} else { cat('Regional inequality ranking table was not available.') }",
        "```",
        "",
        "## 9.4 Lorenz Curve (Geographic Distribution)",
        "```{r lorenz_report_plot, echo=FALSE, fig.width=5, fig.height=5}",
        "if(!is.null(params$outcome_vec)) {",
        "  plot(ineq::Lc(params$outcome_vec), col = '#c0392b', lwd = 3, main = 'Geographic Lorenz Curve')",
        "  abline(0, 1, lty = 2, col = 'gray')",
        "} else { cat('Outcome data not available for Lorenz plot.') }",
        "```",
        "",
        "## 9.5 Policy Implications",
        "A high Gini coefficient (typically > 0.4) indicates that development gains are highly concentrated, suggesting a need for 'Area-Based Targeting' rather than universal policies. If the Indirect Impacts in Section 8.4 are significant, it further suggests that interventions in one region will have positive externalities (spillovers) in neighboring areas.",
        "",
        "\\newpage",
        "# Paper 4: Multilevel Determinants and Variance Decomposition Analysis",
        "## 10.1 Methodology",
        "We employed a 4-stage Hierarchical Linear Model (HLM) to account for the nested structure of the survey data (individuals clustered within communities). The analysis proceeded as follows: (M0) a Null Model to estimate baseline clustering; (M1) an Individual-level model; (M2) a Community-level model; and (M3) a Full Model integrating all covariates. This approach allows us to partition the variance between different levels of the social hierarchy.",
        "",
        "## 10.2 Formal Publication Table: Side-by-Side Model Comparison",
        "The table below presents the fixed effect estimates. For binary or count outcomes, these are reported as Adjusted Odds Ratios (AOR) or Incidence Rate Ratios (ARR) with 95% Confidence Intervals.",
        "```{r ml_table, echo=FALSE}",
        "models <- params$models",
        "tidy_m_full <- function(m, name) {",
        "  if (inherits(m, 'try-error') || is.null(m)) return(data.frame(term='Error', val='-') %>% dplyr::rename(!!name := val))",
        "  res <- broom.mixed::tidy(m, effects = 'fixed', conf.int = TRUE)",
        "  if(params$outcome_type %in% c('Binary', 'Count')) { res <- res %>% mutate(across(c(estimate, conf.low, conf.high), exp)) }",
        "  res %>% mutate(stars = case_when(p.value < 0.001 ~ '***', p.value < 0.01 ~ '**', p.value < 0.05 ~ '*', TRUE ~ '')) %>%",
        "  mutate(formatted = paste0(round(estimate, 2), ' (', round(conf.low, 2), '-', round(conf.high, 2), ')', stars)) %>%",
        "  dplyr::select(term, !!name := formatted)",
        "}",
        "t0 <- tidy_m_full(models$m0, 'M0 (Null)')",
        "t1 <- tidy_m_full(models$m1, 'M1 (Indiv)')",
        "t2 <- tidy_m_full(models$m2, 'M2 (Comm)')",
        "t3 <- tidy_m_full(models$m3, 'M3 (Full)')",
        "final_tab <- t0 %>% full_join(t1, by='term') %>% full_join(t2, by='term') %>% full_join(t3, by='term') %>% filter(!str_detect(term, 'Intercept'))",
        "knitr::kable(final_tab, caption = 'Hierarchical Regression Results (Fixed Effects)')",
        "```",
        "",
        "## 10.3 Model Fitness and Variance Analysis",
        "The Intraclass Correlation Coefficient (ICC) measures the proportion of total variation attributable to community-level factors. The Proportional Change in Variance (PCV) indicates how much of the initial community-level variation was explained by the introduction of covariates.",
        "```{r metrics_table, echo=FALSE}",
        "comp_metrics <- data.frame(",
        "  Model = c('M0 (Null)', 'M1 (Indiv)', 'M2 (Comm)', 'M3 (Full)'),",
        "  AIC = c(AIC(models$m0), AIC(models$m1), AIC(models$m2), AIC(models$m3)),",
        "  BIC = c(BIC(models$m0), BIC(models$m1), BIC(models$m2), BIC(models$m3)),",
        "  ICC = c(calc_icc_internal(models$m0), calc_icc_internal(models$m1), calc_icc_internal(models$m2), calc_icc_internal(models$m3)),",
        "  PCV = c(0, calc_pcv_internal(models$m0, models$m1), calc_pcv_internal(models$m0, models$m2), calc_pcv_internal(models$m0, models$m3))",
        ")",
        "comp_metrics$PCV <- scales::percent(comp_metrics$PCV, accuracy = 0.1)",
        "knitr::kable(comp_metrics, digits = 3, caption = 'Model Selection Criteria and Variance Partitioning')",
        "```",
        "",
        "## 10.4 Visualizing Predictor Strength",
        "The forest plot below visualizes the effect sizes and precision (95% CI) of the predictors in the final model (M3).",
        "```{r forest_plot_rep, echo=FALSE, fig.height=6}",
        "res_f <- broom.mixed::tidy(models$m3, effects = 'fixed', conf.int = TRUE) %>% filter(term != '(Intercept)')",
        "ggplot(res_f, aes(x=estimate, y=reorder(term, estimate))) + ",
        "  geom_vline(xintercept = ifelse(params$outcome_type %in% c('Binary', 'Count'), 1, 0), linetype='dashed', color='red') + ",
        "  geom_errorbarh(aes(xmin=conf.low, xmax=conf.high), height=0.2, color='#2c3e50') + ",
        "  geom_point(size=3, color='#006400') + ",
        "  theme_minimal() + labs(title='Full Model: Fixed Effect Estimates', x='Estimate (Ratio/Beta)', y='')",
        "```",
        "",
        "\\newpage",
        "## 10.5 Spatial Random Effects Map",
        "```{r re_map_rep, echo=FALSE}",
        "re_df <- as.data.frame(ranef(models$m3)$cond$cluster_idx)",
        "re_df$ID <- rownames(re_df)",
        "int_col <- names(re_df)[1]",
        "cluster_region_map <- params$clean_df %>% group_by(!!sym(params$v_cluster)) %>% summarise(Region_Name = first(!!sym(params$v_reg))) %>% mutate(ID = as.character(!!sym(params$v_cluster)))",
        "re_mapped <- re_df %>% left_join(cluster_region_map, by = 'ID') %>% group_by(Region_Name) %>% summarise(mean_re = mean(!!sym(int_col), na.rm=TRUE)) %>% mutate(Match_Key = str_to_upper(str_trim(as.character(Region_Name))))",
        "map_df <- params$spatial_data %>% dplyr::select(Match_Key, geometry) %>% left_join(re_mapped, by = 'Match_Key')",
        "ggplot(map_df) + ",
        "  geom_sf(aes(fill = mean_re), color = 'white', size = 0.1) + ",
        "  scale_fill_gradient2(low='#2980b9', mid='white', high='#c0392b', name='Residual Risk') + ",
        "  theme_void() + labs(title='Map of Unexplained Spatial Variation', subtitle='Aggregated cluster-level random intercepts')",
        "```",
        "",
        "## 10.6 Interpretation",
        "Regions shaded in red indicate areas where the outcome is significantly higher than expected based on the population's characteristics. These areas represent priority zones for further qualitative investigation into unobserved structural barriers.",
        "\\newpage",
        "# 4. Spatiotemporal Dynamics",
        "```{r st_check, include=FALSE}",
        "df_st <- params$st_data",
        "if(is.null(df_st) || nrow(df_st) == 0) { time_points_num <- numeric(0); multi_wave <- FALSE } else { time_points_num <- df_st$time_var %>% as.character() %>% as.numeric() %>% unique() %>% sort(); multi_wave <- length(time_points_num) >= 2 }",
        "shp_geom <- params$spatial_data %>% dplyr::select(Match_Key, geometry) %>% distinct()",
        "```",
        "",
        "\\newpage",
        "# Paper 5: Decomposing the Drivers of Regional Inequality",
        "## 11.1 Methodology",
        "To understand the drivers of geographic inequality, we applied a Blinder-Oaxaca multivariate decomposition analysis. This technique partitions the outcome gap between two geographic groups (e.g., Hot Spots vs. Cold Spots or High-Risk vs. Low-Risk regions) into two distinct components: (1) the **Endowment Effect**, which measures the portion of the gap attributable to differences in population characteristics (e.g., higher education levels in one region); and (2) the **Coefficient Effect**, which measures the portion of the gap due to differences in the 'returns' or effectiveness of those characteristics, often interpreted as structural or environmental advantages.",
        "",
        "## 11.2 Descriptive Profile of Geographic Groups",
        "The table below compares the distribution of key predictors between the two geographic groups. Significant differences in these distributions form the basis of the Endowment Effect.",
        "```{r spatial_t2_rep, echo=FALSE}",
        "if(!is.null(params$spatial_decomp_data)) { ",
        "  knitr::kable(params$spatial_decomp_data$table2, digits=1, caption = 'Table 2: Weighted Characteristics by Region Group') ",
        "} else { cat('Regional decomposition data not available.') }",
        "```",
        "",
        "## 11.3 Outcome Rates by Characteristics",
        "This table shows the actual prevalence of the outcome across different categories within each region. Disparities in these rates across groups contribute to the Coefficient Effect.",
        "```{r spatial_t3_rep, echo=FALSE}",
        "if(!is.null(params$spatial_decomp_data)) { ",
        "  knitr::kable(params$spatial_decomp_data$table3, digits=4, caption = 'Table 3: Outcome Rates by Predictor Category') ",
        "}",
        "```",
        "",
        "## 11.4 Summary Decomposition: Partitioning the Regional Gap",
        "The summary table identifies the 'Bottom Line' of the inequality. It shows the total gap and how much is explained by population characteristics versus structural effects.",
        "```{r spatial_t4_rep, echo=FALSE}",
        "if(!is.null(params$spatial_decomp_data)) { ",
        "  knitr::kable(params$spatial_decomp_data$summary, digits=4, caption = 'Table 4: Summary of Gap Partitioning') ",
        "}",
        "```",
        "",
        "## 11.5 High-Level Drivers of Inequality",
        "The figure below summarizes the percentage contribution of the Explained (Endowment) and Unexplained (Coefficient) components to the total geographic gap.",
        "```{r sp_summary_pct_plot, echo=FALSE, fig.width=6, fig.height=4}",
        "if(!is.null(params$spatial_decomp_data)) {",
        "  df_sum <- params$spatial_decomp_data$summary[4:5,]",
        "  df_sum$Pct <- (df_sum$Value / sum(df_sum$Value)) * 100",
        "  ggplot(df_sum, aes(x = Component, y = Pct, fill = Component)) +",
        "    geom_bar(stat = 'identity', width = 0.6) +",
        "    geom_text(aes(label = paste0(round(Pct, 1), '%')), vjust = -0.5, fontface='bold') +",
        "    scale_fill_manual(values = c('#34495e', '#e67e22')) +",
        "    theme_minimal() + labs(title = 'Total Contribution to Regional Gap (%)', y = 'Percentage Share', x = '') +",
        "    theme(legend.position = 'none')",
        "}",
        "```",
        "",
        "## 11.6 Visualizing Variable-Specific Contributions (Forest Plot)",
        "This forest plot identifies which specific variables are the strongest drivers of inequality. Variables on the right side of the zero line contribute to widening the gap, while those on the left contribute to narrowing it.",
        "```{r sp_fig5_rep, echo=FALSE, fig.width=10, fig.height=7}",
        "if(!is.null(params$spatial_decomp_data)) {",
        "  df_p <- params$spatial_decomp_data$plot_data",
        "  ggplot(df_p, aes(x = Est, y = reorder(Variable, Est))) +",
        "    geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +",
        "    geom_errorbarh(aes(xmin = Est - 1.96*SE, xmax = Est + 1.96*SE), height = 0.2, color = 'gray') +",
        "    geom_point(aes(color = P < 0.05), size = 3) +",
        "    facet_wrap(~ Type, scales = 'free_x') +",
        "    scale_color_manual(values = c('TRUE' = '#c0392b', 'FALSE' = '#2c3e50'), name = 'Statistically Significant') +",
        "    theme_minimal(base_size = 12) + labs(x = 'Contribution to Regional Gap', y = '')",
        "}",
        "```",
        "",
        "## 11.7 Detailed Variable Contributions (Aggregated)",
        "The final table provides the precise numerical contributions of each aggregated variable to the regional gap.",
        "```{r sp_t56_rep, echo=FALSE}",
        "if(!is.null(params$spatial_decomp_data)) { ",
        "  params$spatial_decomp_data$detail %>% ",
        "    dplyr::select(Variable, Endowment, Coefficient, Pct_Share) %>%",
        "    knitr::kable(digits = 4, caption = 'Table 5: Detailed Aggregated Contributions to Inequality') ",
        "}",
        "```",
        "\\newpage",
        "# Paper 6: Spatiotemporal Dynamics and Emerging Hotspots",
        "## 12.1 Methodology",
        "This paper analyzes the evolution of **`r params$v_dv`** across multiple survey waves. It assesses the stability of spatial clustering using Global Moran's I trends, small-multiple maps, rate-of-change maps, LISA transitions, and predictive spatio-temporal interaction modeling.",
        "",
        "## 12.2 Global Spatial Autocorrelation Trend",
        "The figure below tracks the strength of geographic clustering over time. A rising Moran's I suggests that the outcome is becoming more concentrated in specific regions, while a declining trend suggests geographic convergence.",
        "```{r st_moran_trend, echo=FALSE}",
        "if (multi_wave) {",
        "  moran_results <- df_st %>%",
        "    group_by(time_var) %>%",
        "    group_modify(~ {",
        "      current_df <- .x",
        "      nb_st <- poly2nb(current_df, queen = TRUE)",
        "      if(any(card(nb_st) == 0)) nb_st <- knn2nb(knearneigh(st_coordinates(st_centroid(current_df)), k = 1))",
        "      lw_st <- nb2listw(nb_st, style = 'W', zero.policy = TRUE)",
        "      if (length(lw_st$neighbours) > 1) {",
        "        m_test <- moran.test(current_df$outcome, lw_st, zero.policy = TRUE)",
        "        tibble(Moran_I = m_test$estimate[1], P_Value = m_test$p.value)",
        "      } else { tibble(Moran_I = NA, P_Value = NA) }",
        "    }) %>% ungroup() %>% mutate(Significant = P_Value < 0.05)",
        "  ggplot(moran_results, aes(x = time_var, y = Moran_I, group = 1)) +",
        "    geom_line(color = '#1e8449', size = 1.5) +",
        "    geom_point(aes(color = Significant), size = 4) +",
        "    geom_hline(yintercept = 0, linetype = 'dashed', color = 'gray') +",
        "    scale_color_manual(values = c('TRUE' = '#c0392b', 'FALSE' = '#bdc3c7')) +",
        "    labs(title = 'Evolution of Spatial Clustering (Global Moran\\'s I)', x = 'Survey Wave', y = 'Moran\\'s I') +",
        "    theme_minimal()",
        "} else { cat('Skipped: Requires multiple time points.') }",
        "```",
        "",
        "## 12.3 Spatiotemporal Evolution (Small Multiples)",
        "The maps below visualize the regional distribution of the outcome across all available survey waves.",
        "```{r st_small_multiples, echo=FALSE, fig.height=8}",
        "if (!is.null(df_st) && nrow(df_st) > 0) {",
        "  ggplot(df_st) +",
        "    geom_sf(aes(fill = outcome), color = 'white', size = 0.05) +",
        "    scale_fill_viridis_c(option = 'plasma', labels = scales::percent) +",
        "    facet_wrap(~ time_var, ncol = 3) +",
        "    labs(title = 'Regional Prevalence Rates by Time Period', fill = 'Rate') +",
        "    theme_void() + theme(strip.text = element_text(face = 'bold', size = 11))",
        "}",
        "```",
        "",
        "## 12.4 Analysis of Change: Rate of Change and LISA Transitions",
        "The 'Rate of Change' map identifies regions of improvement (blue) versus decline (red). The 'LISA Transition' map tracks the movement of hotspots (HH) and coldspots (LL) between the first and last wave.",
        "```{r st_dynamics_combined, echo=FALSE, fig.width=10, fig.height=5}",
        "if (multi_wave) {",
        "  t_min <- min(time_points_num); t_max <- max(time_points_num)",
        "  df_min <- df_st %>% filter(as.numeric(as.character(time_var)) == t_min) %>% st_drop_geometry() %>% dplyr::select(Match_Key, outcome_min = outcome)",
        "  df_max <- df_st %>% filter(as.numeric(as.character(time_var)) == t_max) %>% st_drop_geometry() %>% dplyr::select(Match_Key, outcome_max = outcome)",
        "  df_change <- df_min %>% full_join(df_max, by = 'Match_Key') %>% mutate(pct_change = (outcome_max - outcome_min) / outcome_min * 100)",
        "  map_change <- shp_geom %>% left_join(df_change, by = 'Match_Key')",
        "  p_roc <- ggplot(map_change) + geom_sf(aes(fill = pct_change), color = 'white', size = 0.1) +",
        "    scale_fill_gradient2(low = '#2980b9', mid = 'white', high = '#c0392b', name = '% Change') +",
        "    theme_void() + labs(title = 'Rate of Change')",
        "  print(p_roc)",
        "}",
        "```",
        "",
        "## 12.5 APA Table: Moran's I Trend and P-values",
        "```{r st_moran_table_rep, echo=FALSE}",
        "if(!is.null(params$st_moran_table)) { knitr::kable(params$st_moran_table, caption = 'APA Table: Spatiotemporal Moran Trend by Survey Wave') } else { cat('Moran trend table not available.') }",
        "```",
        "",
        "## 12.6 APA Table: Rate of Change by Region",
        "```{r st_rate_table_rep, echo=FALSE}",
        "if(!is.null(params$st_rate_table)) { knitr::kable(params$st_rate_table, digits = 4, caption = 'APA Table: Rate of Change by Region') } else { cat('Rate of change table not available.') }",
        "```",
        "",
        "## 12.7 APA Table: LISA Transition Classification",
        "```{r st_lisa_table_rep, echo=FALSE}",
        "if(!is.null(params$st_lisa_table)) { knitr::kable(params$st_lisa_table, caption = 'APA Table: LISA Transition Classification') } else { cat('LISA transition table not available.') }",
        "```",
        "",
        "## 12.8 Predictive Spatio-Temporal Interaction Model",
        "Beyond descriptive trends, this model partitions geographic risk into a global temporal trend and region-specific deviations. It identifies emerging hotspots where the local trend worsens faster than the overall trend.",
        "",
        "### 12.8.1 Interaction Model Results with P-values",
        "```{r st_int_coef_rep_paper6, echo=FALSE}",
        "if(!is.null(params$st_int_model)) {",
        "  res_int <- tidy_glmmTMB_fixed_apa(params$st_int_model, exponentiate = TRUE)",
        "  knitr::kable(res_int, digits = 3, caption = 'APA Table: Spatiotemporal Interaction Model Fixed Effects')",
        "} else { cat('Interaction model was not run.') }",
        "```",
        "",
        "### 12.8.2 Emerging Hotspot Ranking",
        "```{r st_hotspot_table_rep, echo=FALSE}",
        "if(!is.null(params$st_hotspot_table)) { knitr::kable(params$st_hotspot_table, digits = 4, caption = 'APA Table: Emerging Hotspot Ranking') } else { cat('Emerging hotspot table not available.') }",
        "```",
        "",
        "### 12.8.3 Mapping Emerging Hotspots",
        "```{r st_int_map_rep_paper6, echo=FALSE}",
        "if(!is.null(params$st_int_data)) {",
        "  map_int <- shp_geom %>% left_join(params$st_int_data, by = 'Match_Key')",
        "  ggplot(map_int) + geom_sf(aes(fill = interaction_effect), color = 'white', size = 0.1) + scale_fill_gradient2(low = '#2980b9', mid = 'white', high = '#c0392b', name = 'Interaction') + theme_void() + labs(title = 'Map of Emerging Hotspots', subtitle = 'Positive values indicate rising relative risk')",
        "} else { cat('Interaction data not available for mapping.') }",
        "```",
        "",
        "\\newpage",
        "# Paper 7: Spatiotemporal Decomposition of Change",
        "## 13.1 Methodology",
        "This paper applies spatial and temporal Blinder-Oaxaca decomposition to partition the change between the baseline and endline survey waves into Explained (Endowment) and Unexplained (Coefficient) components.",
        "",
        "## 13.2 Summary Decomposition: Partitioning the Time Gap",
        "This table summarizes the total change over time and partitions it into the portion explained by shifts in population characteristics (Endowment) versus shifts in the effectiveness of those characteristics (Coefficient).",
        "```{r st_t4_rep, echo=FALSE}",
        "if(!is.null(params$decomp_data)) { ",
        "  knitr::kable(params$decomp_data$summary, digits=4, caption = 'Table 4: Summary of Spatiotemporal Gap Partitioning') ",
        "}",
        "```",
        "",
        "## 13.3 High-Level Drivers of Change (%)",
        "```{r st_summary_pct_plot, echo=FALSE, fig.width=6, fig.height=4}",
        "if(!is.null(params$decomp_data)) {",
        "  df_sum <- params$decomp_data$summary[4:5,]",
        "  df_sum$Pct <- (df_sum$Value / sum(df_sum$Value)) * 100",
        "  ggplot(df_sum, aes(x = Component, y = Pct, fill = Component)) +",
        "    geom_bar(stat = 'identity', width = 0.6) +",
        "    geom_text(aes(label = paste0(round(Pct, 1), '%')), vjust = -0.5, fontface='bold') +",
        "    scale_fill_manual(values = c('#2c3e50', '#27ae60')) +",
        "    theme_minimal() + labs(title = 'Total Contribution to Change (%)', y = 'Percentage Share', x = '') +",
        "    theme(legend.position = 'none')",
        "}",
        "```",
        "",
        "## 13.4 Visualizing Drivers of Change (Forest Plot)",
        "This plot identifies which specific variables were the strongest drivers of the observed change between waves.",
        "```{r st_fig5_rep, echo=FALSE, fig.width=10, fig.height=7}",
        "if(!is.null(params$decomp_data)) {",
        "  df_p <- params$decomp_data$plot_data",
        "  ggplot(df_p, aes(x = Est, y = reorder(Variable, Est))) +",
        "    geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +",
        "    geom_errorbarh(aes(xmin = Est - 1.96*SE, xmax = Est + 1.96*SE), height = 0.2, color = 'gray') +",
        "    geom_point(aes(color = P < 0.05), size = 3) +",
        "    facet_wrap(~ Type, scales = 'free_x') +",
        "    scale_color_manual(values = c('TRUE' = '#c0392b', 'FALSE' = '#2c3e50'), name = 'Significant') +",
        "    theme_minimal(base_size = 12) + labs(x = 'Contribution to Change', y = '')",
        "}",
        "```",
        "",
        "## 13.5 Detailed Variable Contributions (Aggregated)",
        "The table below provides the precise numerical contributions of each aggregated variable to the change over time.",
        "```{r st_t56_rep, echo=FALSE}",
        "if(!is.null(params$decomp_data)) { ",
        "  params$decomp_data$detail %>% ",
        "    dplyr::select(Variable, Endowment, Coefficient, Pct_Share) %>%",
        "    knitr::kable(digits = 4, caption = 'Table 5: Detailed Aggregated Contributions to Change') ",
        "}",
        "```",
        "",
        "\\newpage",
        "# Paper 8: Joint Spatial Modeling and Shared Determinants",
        "## 13.1 Methodology",
        "To investigate the interconnectedness of the two selected outcomes, we employed a Frequentist Joint Multivariate Generalized Linear Mixed Model (GLMM) using the `glmmTMB` engine. This approach allows for the simultaneous estimation of determinants for both outcomes while accounting for shared unexplained variance at the cluster level. By modeling the outcomes jointly, we can identify whether specific socio-economic or environmental factors act as shared drivers for multiple public health vulnerabilities.",
        "",
        "## 13.2 Shared Variance and Spatial Clustering",
        "The table below reports the cluster-level correlation between the two outcomes and the Global Moran's I for the joint burden (the geographic overlap where both conditions are present).",
        "```{r joint_metrics_rep, echo=FALSE}",
        "if(!is.null(params$joint_model)) {",
        "  # 1. Extract Shared Correlation",
        "  vc <- VarCorr(params$joint_model)",
        "  corr_val <- attr(vc$cond$cluster_idx, 'correlation')[1,2]",
        "  ",
        "  # 2. Calculate Joint Moran's I",
        "  df_js <- params$joint_spatial",
        "  nb_js <- poly2nb(df_js, queen = TRUE)",
        "  if(any(card(nb_js) == 0)) nb_js <- knn2nb(knearneigh(st_coordinates(st_centroid(df_js)), k = 1))",
        "  lw_js <- nb2listw(nb_js, style = 'W', zero.policy = TRUE)",
        "  m_js <- moran.test(df_js$joint_rate, lw_js, zero.policy = TRUE)",
        "  ",
        "  joint_metrics <- data.frame(",
        "    Metric = c('Shared Cluster Correlation (Rho)', 'Global Moran\\'s I (Joint Burden)', 'P-value (Spatial Clustering)'),",
        "    Value = c(round(corr_val, 3), round(m_js$estimate[1], 3), format.pval(m_js$p.value))",
        "  )",
        "  knitr::kable(joint_metrics, caption = 'Shared Variance and Spatial Autocorrelation of Joint Burden')",
        "} else { cat('Joint modeling was not performed.') }",
        "```",
        "",
        "## 13.3 Joint Determinants: Shared Drivers of Vulnerability",
        "The forest plot visualizes the effect sizes (Log-Odds) for predictors across both outcomes. Predictors that show significant effects in the same direction for both outcomes represent shared structural drivers.",
        "```{r joint_forest_rep, echo=FALSE, fig.width=10, fig.height=7}",
        "if(!is.null(params$joint_model)) {",
        "  res_j <- broom.mixed::tidy(params$joint_model, effects = 'fixed', conf.int = TRUE) %>%",
        "    filter(!term %in% c('outcome_typeoutcome1', 'outcome_typeoutcome2')) %>%",
        "    mutate(",
        "      Outcome = ifelse(str_detect(term, 'outcome1'), 'Outcome 1', 'Outcome 2'),",
        "      Variable = str_replace(term, 'outcome_typeoutcome[12]:', '')",
        "    )",
        "  ",
        "  ggplot(res_j, aes(x = estimate, y = Variable, color = Outcome)) +",
        "    geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +",
        "    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.3, position = position_dodge(width = 0.5)) +",
        "    geom_point(size = 3, position = position_dodge(width = 0.5)) +",
        "    scale_color_manual(values = c('Outcome 1' = '#2c3e50', 'Outcome 2' = '#e67e22')) +",
        "    theme_minimal(base_size = 12) +",
        "    labs(title = 'Joint Model Forest Plot', subtitle = 'Comparing predictor strengths across both outcomes', x = 'Log-Odds Ratio', y = '')",
        "}",
        "```",
        "",
        "## 13.4 Joint Hotspot Analysis (LISA Map)",
        "The map below identifies 'Joint Hotspots' (High-High clusters), representing geographic areas where both outcomes are significantly higher than the national average and are surrounded by similarly high-burden areas.",
        "```{r joint_lisa_rep, echo=FALSE}",
        "if(!is.null(params$joint_spatial)) {",
        "  df_l <- params$joint_spatial",
        "  nb_l <- poly2nb(df_l, queen = TRUE)",
        "  if(any(card(nb_l) == 0)) nb_l <- knn2nb(knearneigh(st_coordinates(st_centroid(df_l)), k = 1))",
        "  lw_l <- nb2listw(nb_l, style = 'W', zero.policy = TRUE)",
        "  loc_l <- localmoran(df_l$joint_rate, lw_l, zero.policy = TRUE)",
        "  m_avg_l <- mean(df_l$joint_rate, na.rm=TRUE)",
        "  df_l$Pr <- loc_l[,5]",
        "  df_l$lag_j <- lag.listw(lw_l, df_l$joint_rate, zero.policy = TRUE)",
        "  df_l$lisa <- case_when(",
        "    df_l$Pr <= 0.05 & df_l$joint_rate >= m_avg_l & df_l$lag_j >= m_avg_l ~ 'High-High (Hotspot)',",
        "    df_l$Pr <= 0.05 & df_l$joint_rate <= m_avg_l & df_l$lag_j <= m_avg_l ~ 'Low-Low (Coldspot)',",
        "    TRUE ~ 'Not Significant'",
        "  )",
        "  ggplot(df_l) + geom_sf(aes(fill = lisa), color = 'white', size = 0.1) +",
        "    scale_fill_manual(values = c('High-High (Hotspot)'='#c0392b', 'Low-Low (Coldspot)'='#2980b9', 'Not Significant'='#ecf0f1')) +",
        "    theme_void() + labs(title = 'Joint LISA Cluster Map', subtitle = 'Geographic overlap of high-burden areas', fill = 'Cluster Type')",
        "}",
        "```",
        "",
        "## 13.5 Joint Model Coefficients Table (Adjusted Odds Ratios)",
        "The table below provides the formal regression results. Estimates are reported as Adjusted Odds Ratios (AOR).",
        "```{r joint_table_rep, echo=FALSE}",
        "if(!is.null(params$joint_model)) {",
        "  res_tab <- broom.mixed::tidy(params$joint_model, effects = 'fixed', conf.int = TRUE) %>%",
        "    mutate(across(c(estimate, conf.low, conf.high), exp)) %>%",
        "    mutate(P_Value = format.pval(p.value, digits = 3)) %>%",
        "    dplyr::select(Term = term, AOR = estimate, LCI = conf.low, UCI = conf.high, P_Value)",
        "  knitr::kable(res_tab, digits = 3, caption = 'Multivariate Joint Model Coefficients (AOR)')",
        "}",
        "```",
        "\\newpage",
        "# Paper 9: Decomposing the Regional Gap in the Syndemic Burden",
        "## 14.1 Methodology",
        "This paper applies a Blinder-Oaxaca multivariate decomposition to the **joint occurrence** of Outcome 1 and Outcome 2. We analyze why the 'Syndemic Burden' (children suffering from both conditions) differs between geographic regions.",
        "",
        "## 14.2 Table 2: Descriptive Profile of Joint Groups",
        "```{r jd_t2_rep, echo=FALSE}",
        "if(!is.null(params$joint_decomp_data)) {",
        "  knitr::kable(params$joint_decomp_data$table2, digits=1, caption = 'Weighted Characteristics by Region Group')",
        "}",
        "```",
        "",
        "## 14.3 Table 3: Joint Outcome Rates by Characteristics",
        "```{r jd_t3_rep, echo=FALSE}",
        "if(!is.null(params$joint_decomp_data)) {",
        "  knitr::kable(params$joint_decomp_data$table3, digits=4, caption = 'Joint Outcome Rates by Predictor Category')",
        "}",
        "```",
        "",
        "## 14.4 Table 4: Summary of Joint Gap Partitioning",
        "```{r jd_t4_rep, echo=FALSE}",
        "if(!is.null(params$joint_decomp_data)) {",
        "  knitr::kable(params$joint_decomp_data$summary, digits=4, caption = 'Partitioning the Regional Gap in Joint Burden')",
        "}",
        "```",
        "",
        "## 14.5 Figure 5: Drivers of the Joint Gap (Forest Plot)",
        "```{r jd_fig5_rep, echo=FALSE, fig.width=10, fig.height=7}",
        "if(!is.null(params$joint_decomp_data)) {",
        "  df_p8 <- params$joint_decomp_data$detail",
        "  ggplot(df_p8, aes(x = Endowment, y = reorder(Variable, Endowment))) +",
        "    geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +",
        "    geom_point(size = 3, color = '#c0392b') +",
        "    theme_minimal(base_size = 12) +",
        "    labs(title = 'Variable Contributions to the Joint Burden Gap', x = 'Contribution Size', y = '')",
        "}",
        "```",
        "",
        "## 14.6 Table 5: Detailed Aggregated Contributions",
        "```{r jd_t5_rep, echo=FALSE}",
        "if(!is.null(params$joint_decomp_data)) {",
        "  knitr::kable(params$joint_decomp_data$detail, digits=4, caption = 'Aggregated Variable Contributions to the Joint Gap')",
        "}",
        "```",
        "",
        "## 14.7 Policy Conclusion",
        "If the **Coefficient Effect** is dominant, it suggests that even if population characteristics were equalized, the regional gap would persist due to structural inequalities (e.g., climate, health infrastructure). If the **Endowment Effect** is dominant, policies should focus on improving individual-level factors like education and wealth in the high-burden region.",
        "\\newpage",
        "# Paper 10: Spatiotemporal Drivers of the Syndemic Burden",
        "## 15.1 Methodology",
        "This paper analyzes the **evolution** of the joint burden of Outcome 1 and Outcome 2 over time. We decompose the change between Wave 1 and Wave 2 into population shifts (Endowment) and structural shifts (Coefficient).",
        "",
        "## 15.2 Table 2: Descriptive Profile of Waves",
        "```{r jsd_t2_rep, echo=FALSE}",
        "if(!is.null(params$joint_st_decomp_data)) {",
        "  knitr::kable(params$joint_st_decomp_data$table2, digits=1, caption = 'Weighted Characteristics by Wave')",
        "}",
        "```",
        "",
        "## 15.3 Table 4: Summary of Joint Change Partitioning",
        "```{r jsd_t4_rep, echo=FALSE}",
        "if(!is.null(params$joint_st_decomp_data)) {",
        "  knitr::kable(params$joint_st_decomp_data$summary, digits=4, caption = 'Partitioning the Change in Joint Burden Over Time')",
        "}",
        "```",
        "",
        "## 15.4 Figure 5: Drivers of Joint Change (Forest Plot)",
        "```{r jsd_fig5_rep, echo=FALSE, fig.width=10, fig.height=7}",
        "if(!is.null(params$joint_st_decomp_data)) {",
        "  df_p9 <- params$joint_st_decomp_data$detail",
        "  ggplot(df_p9, aes(x = Endowment, y = reorder(Variable, Endowment))) +",
        "    geom_vline(xintercept = 0, linetype = 'dashed', color = 'red') +",
        "    geom_point(size = 3, color = '#27ae60') +",
        "    theme_minimal(base_size = 12) +",
        "    labs(title = 'Variable Contributions to the Change in Joint Burden', x = 'Contribution Size', y = '')",
        "}",
        "```",
        "\\newpage",
        "# Paper 11: Spatial Accessibility and Service Desert Analysis",
        "## 16.1 Methodology",
        "This paper examines whether geographic access to services is associated with spatial differences in the outcome. Where service-point data are available, the module computes distance-based access measures. Where service-point data are not uploaded, it uses the selected accessibility proxy as a regional access indicator.",
        "",
        "## 16.2 APA Accessibility Summary by Region",
        "```{r p11_t1_report, echo=FALSE}",
        "if(!is.null(params$paper11_data)) { knitr::kable(params$paper11_data$summary, digits = 4, caption = 'APA Table: Accessibility Summary by Region') } else { cat('Paper 11 analysis was not run.') }",
        "```",
        "",
        "## 16.3 Service Desert Classification",
        "```{r p11_t2_report, echo=FALSE}",
        "if(!is.null(params$paper11_data)) { knitr::kable(params$paper11_data$desert, digits = 4, caption = 'APA Table: Service Desert Classification') }",
        "```",
        "",
        "## 16.4 Outcome Model by Accessibility",
        "```{r p11_t3_report, echo=FALSE}",
        "if(!is.null(params$paper11_data)) { knitr::kable(params$paper11_data$model, digits = 4, caption = 'APA Table: Outcome Model by Accessibility') }",
        "```",
        "",
        "## 16.5 Figures: Service Desert Map and Access-Outcome Gradient",
        "```{r p11_figs_report, echo=FALSE, fig.width=10, fig.height=5}",
        "if(!is.null(params$paper11_data)) {",
        "  p1 <- ggplot(params$paper11_data$data) + geom_sf(aes(fill = Desert_Status), color = 'white', size = 0.1) + scale_fill_manual(values = c('Service desert / low access' = '#b2182b', 'Moderate / high access' = '#2166ac'), na.value = 'grey80') + theme_void() + labs(title = 'Service Desert Classification', fill = 'Status')",
        "  d11 <- sf::st_drop_geometry(params$paper11_data$data)",
        "  p2 <- ggplot(d11, aes(x = access_metric, y = outcome)) + geom_point(size = 3, color = '#006400') + geom_smooth(method = 'lm', se = TRUE, color = '#c0392b') + theme_minimal() + labs(title = 'Accessibility-Outcome Gradient', x = 'Access metric', y = 'Outcome')",
        "  print(cowplot::plot_grid(p1, p2, ncol = 2))",
        "}",
        "```",
        "",
        "\\newpage",
        "# Paper 12: Geographically Weighted / Local Determinants Analysis",
        "## 17.1 Methodology",
        "This paper assesses whether predictor-outcome relationships vary across space. The module estimates global and local neighborhood-specific slopes and maps local coefficients and local model fit.",
        "",
        "## 17.2 Global vs Local Model Summary",
        "```{r p12_t1_report, echo=FALSE}",
        "if(!is.null(params$paper12_data)) { knitr::kable(params$paper12_data$global, digits = 4, caption = 'APA Table: Global and Local Model Summary') } else { cat('Paper 12 analysis was not run.') }",
        "```",
        "",
        "## 17.3 Local Coefficients by Region",
        "```{r p12_t2_report, echo=FALSE}",
        "if(!is.null(params$paper12_data)) { knitr::kable(params$paper12_data$local, digits = 4, caption = 'APA Table: Local Coefficients by Region') }",
        "```",
        "",
        "## 17.4 Figures: Local Coefficient and Local R-squared Maps",
        "```{r p12_figs_report, echo=FALSE, fig.width=10, fig.height=5}",
        "if(!is.null(params$paper12_data)) {",
        "  p1 <- ggplot(params$paper12_data$data) + geom_sf(aes(fill = Local_Coefficient), color = 'white', size = 0.1) + scale_fill_gradient2(low = '#2980b9', mid = 'white', high = '#c0392b') + theme_void() + labs(title = 'Local Coefficient Map', fill = 'Local beta')",
        "  p2 <- ggplot(params$paper12_data$data) + geom_sf(aes(fill = Local_R2), color = 'white', size = 0.1) + scale_fill_viridis_c(option = 'mako') + theme_void() + labs(title = 'Local Model Fit Map', fill = 'Local R-squared')",
        "  print(cowplot::plot_grid(p1, p2, ncol = 2))",
        "}",
        "```",
        "",
        "\\newpage",
        "# Paper 13: Spatial Priority Index and Targeting Analysis",
        "## 18.1 Methodology",
        "This paper develops a composite spatial priority index by combining the mapped outcome with selected vulnerability indicators. The purpose is to rank areas for targeted planning and resource allocation.",
        "",
        "## 18.2 Spatial Priority Ranking",
        "```{r p13_t1_report, echo=FALSE}",
        "if(!is.null(params$paper13_data)) { knitr::kable(params$paper13_data$priority, digits = 4, caption = 'APA Table: Spatial Priority Ranking') } else { cat('Paper 13 analysis was not run.') }",
        "```",
        "",
        "## 18.3 Indicator Weights and Contributions",
        "```{r p13_t2_report, echo=FALSE}",
        "if(!is.null(params$paper13_data)) { knitr::kable(params$paper13_data$indicators, digits = 4, caption = 'APA Table: Indicator Contribution Summary') }",
        "```",
        "",
        "## 18.4 Figures: Priority Index Map and Indicator Contribution Plot",
        "```{r p13_figs_report, echo=FALSE, fig.width=10, fig.height=5}",
        "if(!is.null(params$paper13_data)) {",
        "  p1 <- ggplot(params$paper13_data$data) + geom_sf(aes(fill = Priority_Score), color = 'white', size = 0.1) + scale_fill_viridis_c(option = 'magma', direction = -1) + theme_void() + labs(title = 'Spatial Priority Index', fill = 'Score')",
        "  p2 <- ggplot(params$paper13_data$indicators, aes(x = reorder(Indicator, Mean_Scaled_Score), y = Mean_Scaled_Score)) + geom_col(fill = '#006400') + coord_flip() + theme_minimal() + labs(title = 'Indicator Contribution Summary', x = '', y = 'Mean scaled score')",
        "  print(cowplot::plot_grid(p1, p2, ncol = 2))",
        "}",
        "```",
        "",
        "\\newpage",
        "# Paper 14: Spatial Scan and Cluster Detection Analysis",
        "## 19.1 Methodology",
        "This paper detects statistically meaningful high-burden and cold-cluster areas by combining global autocorrelation, local cluster classification, Getis-Ord logic, and high-burden thresholding.",
        "",
        "## 19.2 Global Cluster Diagnostics",
        "```{r p14_t1_report, echo=FALSE}",
        "if(!is.null(params$paper14_data)) { knitr::kable(params$paper14_data$global, digits = 4, caption = 'APA Table: Global Cluster Diagnostics') } else { cat('Paper 14 analysis was not run.') }",
        "```",
        "",
        "## 19.3 Local Cluster Classification",
        "```{r p14_t2_report, echo=FALSE}",
        "if(!is.null(params$paper14_data)) { knitr::kable(params$paper14_data$local, digits = 4, caption = 'APA Table: Local Cluster Classification') }",
        "```",
        "",
        "## 19.4 Figures: Candidate Spatial Cluster Map and Burden Profile",
        "```{r p14_figs_report, echo=FALSE, fig.width=10, fig.height=5}",
        "if(!is.null(params$paper14_data)) {",
        "  p1 <- ggplot(params$paper14_data$data) + geom_sf(aes(fill = Cluster_Class), color = 'white', size = 0.1) + scale_fill_manual(values = c('High-burden hot cluster' = '#b2182b', 'High-burden candidate' = '#ef8a62', 'Cold cluster' = '#2166ac', 'Not clustered' = '#d9d9d9'), na.value = 'grey80') + theme_void() + labs(title = 'Spatial Cluster Detection Map', fill = 'Cluster')",
        "  p2 <- ggplot(params$paper14_data$profile, aes(x = reorder(Cluster_Class, Mean_Outcome), y = Mean_Outcome)) + geom_col(fill = '#006400') + coord_flip() + theme_minimal() + labs(title = 'Cluster Burden Profile', x = '', y = 'Mean outcome')",
        "  print(cowplot::plot_grid(p1, p2, ncol = 2))",
        "}",
        "```",
        "",
        "\\newpage",
        "# Paper 15: Spatial Equity and Environmental Justice Analysis",
        "## 20.1 Methodology",
        "This paper examines whether vulnerable groups experience higher outcome burdens and whether those burdens are spatially concentrated. It supports equity-focused and justice-oriented targeting.",
        "",
        "## 20.2 Equity Gap Summary",
        "```{r p15_t1_report, echo=FALSE}",
        "if(!is.null(params$paper15_data)) { knitr::kable(params$paper15_data$gap, digits = 4, caption = 'APA Table: Equity Gap Summary') } else { cat('Paper 15 analysis was not run.') }",
        "```",
        "",
        "## 20.3 Regional Equity Burden Table",
        "```{r p15_t2_report, echo=FALSE}",
        "if(!is.null(params$paper15_data)) { knitr::kable(params$paper15_data$region, digits = 4, caption = 'APA Table: Regional Equity Burden') }",
        "```",
        "",
        "## 20.4 Figures: Vulnerable Group Share and Outcome by Equity Group",
        "```{r p15_figs_report, echo=FALSE, fig.width=10, fig.height=5}",
        "if(!is.null(params$paper15_data)) {",
        "  p1 <- ggplot(params$paper15_data$map) + geom_sf(aes(fill = Vulnerable_Share), color = 'white', size = 0.1) + scale_fill_viridis_c(option = 'magma', labels = scales::percent) + theme_void() + labs(title = 'Vulnerable Group Share by Region', fill = 'Share')",
        "  p2 <- ggplot(params$paper15_data$group, aes(x = reorder(eq_group, Outcome), y = Outcome)) + geom_col(fill = '#006400') + coord_flip() + theme_minimal() + labs(title = 'Outcome by Equity Group', x = '', y = 'Weighted outcome')",
        "  print(cowplot::plot_grid(p1, p2, ncol = 2))",
        "}",
        "```",
        "",
        "\\newpage",
        "# Citation and Acknowledgments",
        "If you use **HawaSpatial Pro** in your research, please cite it as:",
        "",
        "Muse, A. H. (2026). HawaSpatial: Holistic and Areal Weighted Analysis for Global Development (Version 0.1.10). Zenodo. https://doi.org/10.5281/zenodo.18290101",
        "",
        "### BibTeX for LaTeX users:",
        "```bibtex",
        "@manual{muse2026hawa,",
        "  title = {HawaSpatial: Holistic and Areal Weighted Analysis for Global Development},",
        "  author = {Muse, Abdisalam Hassan},",
        "  year = {2026},",
        "  note = {R package version 0.1.10},",
        "  doi = {10.5281/zenodo.18290101},",
        "  url = {https://doi.org/10.5281/zenodo.18290101}",
        "}",
        "```",
        ""
      )

      writeLines(report_content, tempReport)

      rmarkdown::render(tempReport, output_file = file,
                        params = list(
                          doc_title = input$rep_title,
                          spatial_data = safe_get(adal_spatial()),
                          models = safe_get(freq_res()),
                          outcome_type = input$outcome_type,
                          v_reg = input$v_reg,
                          v_cluster = input$v_cluster,
                          clean_df = safe_get(processed_df()),
                          raw_df = safe_get(raw_data()),
                          v_indiv_cont = input$v_indiv_cont,
                          v_indiv_fact = input$v_indiv_fact,
                          v_comm_cont = input$v_comm_cont,
                          v_comm_fact = input$v_comm_fact,
                          v_wgt = input$v_wgt,
                          is_dhs_wgt = input$is_dhs_wgt,
                          st_data = safe_get(adal_spatial and temporal()),
                          v_time = input$v_time,
                          v_wealth = input$v_wealth,
                          cov_vars = c(input$v_indiv_cont, input$v_indiv_fact, input$v_comm_cont, input$v_comm_fact),
                          sae_data = safe_get(sae_res()),
                          econ_data = safe_get(econ_res()),
                          econ_type = input$econ_model_type,
                          econ_lm = if(!is.null(safe_get(econ_res()))) econ_res()$lm_tests else NULL,
                          econ_fit_table = safe_get(econ_fit_table_data()),
                          econ_weights_table = safe_get(econ_weights_table_data()),
                          econ_resid_table = safe_get(econ_resid_table_data()),
                          gini_val = if(!is.null(safe_get(ineq_stats()))) ineq_stats()$gini else NULL,
                          conc_val = safe_get(conc_index()),
                          theil_val = if(!is.null(safe_get(ineq_stats()))) ineq_stats()$theil else NULL,
                          outcome_vec = if(!is.null(safe_get(adal_spatial()))) adal_spatial()$outcome else NULL,
                          decomp_data = safe_get(decomp_res()),
                          spatial_decomp_data = safe_get(spatial_decomp_res()),
                          joint_model = safe_get(joint_res()),
                          joint_spatial = safe_get(joint_spatial_data()),
                          joint_decomp_data = safe_get(jd_res()),
                          joint_st_decomp_data = safe_get(jsd_res()),
                          st_int_model = safe_get(st_int_res()$model),
                          st_int_data = safe_get(st_int_res()$interaction_data),
                          st_moran_table = safe_get(st_moran_trend_results()),
                          st_rate_table = safe_get(st_rate_change_table_data()),
                          st_lisa_table = safe_get(st_lisa_transition_table_data()),
                          st_hotspot_table = safe_get(st_hotspot_table_data()),
                          bivariate_table = safe_get(bivariate_screening_table()),
                          ineq_apa_table = if(!is.null(safe_get(inequality_apa_tables()))) inequality_apa_tables()$summary else NULL,
                          ineq_region_table = if(!is.null(safe_get(inequality_apa_tables()))) inequality_apa_tables()$regional else NULL,
                          paper11_data = safe_get(paper11_res()),
                          paper12_data = safe_get(paper12_res()),
                          paper13_data = safe_get(paper13_res()),
                          paper14_data = safe_get(paper14_res()),
                          paper15_data = safe_get(paper15_res())
                        ),
                        envir = new.env(parent = globalenv()))
    }
  )
}


shinyApp(ui, server)
