# =============================================================================
# 05_economic_extension.R
# Group 10 — SMB-A Capstone: Conflict Spillover Analysis
# Team: Karen David, Candace Offman, Mark Oleson
# Advisor: Dr. Elliot Bendoly | OSU Fisher College of Business
# =============================================================================
#
# PURPOSE:
#   Section 5.2b of the S-A3 analysis — extends the NB GLM from 02_predictive.R
#   with three external economic predictors merged at the country-year level:
#
#     (1) World Bank WDI — GDP per capita (constant 2015 USD), GDP growth rate
#         Source: wbstats R package | Indicators: NY.GDP.PCAP.KD, NY.GDP.MKTP.KD.ZG
#
#     (2) Medina & Schneider (2019) Shadow Economy — % of official GDP
#         Source: CESifo WP 7981 / IMF WP/18/17 / Schneider JKU page
#         Adjusts GDP per capita to account for informal economy size.
#         Shadow economy pct also included as an independent predictor
#         (higher shadow economy → weaker state capacity → potential conflict risk).
#
#     (3) IMF Primary Commodity Prices — country-specific commodity price shock
#         Source: IMF PCPS (https://www.imf.org/external/np/res/commod/)
#         Each SCAD country is mapped to its primary commodity export.
#         Shock = annual % change in primary commodity price (year-over-year).
#         Both GDP growth and commodity shock are lagged one year (t-1 → t)
#         to reduce endogeneity (economic conditions predict future conflict,
#         not same-year conflict that may depress economic activity).
#
# MANAGERIAL OBJECTIVE ADDRESSED:
#   MO3: Predict Conflict Severity via Economic Impact
#
# DEPENDS ON:
#   00_data_load_clean.R — combined, model_df (must be sourced first)
#   02_predictive.R      — nb_model, train_df, test_df (for baseline comparison)
#
# OUTPUTS:
#   scad_econ            — combined dataset with economic columns added
#   model_data_econ      — analysis-ready subset (non-missing economic vars)
#   nb_model_econ        — extended NB GLM object
#   scad_wb_crosswalk    — country name crosswalk data frame (also written to CSV)
#   Table E1             — Country crosswalk + missingness report
#   Table E2             — Extended NB GLM coefficient comparison (baseline vs. extended)
#   Table E3             — GVIF table for extended model
#   Table E4             — Model performance comparison (pseudo-R², AIC, N)
#   Figure E1            — Economic predictor coefficient plot (extended model only)
#   Figure E2            — Updated what-if counterfactual (escalation cap at Level 3)
#
# NARRATIVE CAVEATS (appear as Rmd prose, not cat() output):
#   - shadow_pct Route C proxy lacks country-level variation; do not interpret
#     as causal when regional averages are used.
#   - Counterfactual estimates are in-model, not causal guarantees.
#   - Economic model uses a subset of events (missing GDP data excluded);
#     percentage reductions apply to that subset only.
#   - IMF commodity CSV must be downloaded manually if unavailable; see Section 4.
#   - South Sudan: World Bank has no data before 2011 (statehood).
#   - Somalia, Eritrea, Libya: expect substantial gaps in GDP coverage.
#
# DATA PROVENANCE (for Rmd Section 5.2b prose):
#   World Bank World Development Indicators, via wbstats R package.
#   Medina, L. & Schneider, F. (2019). Shedding Light on the Shadow Economy:
#     A Global Database. CESifo WP 7981. SSRN: 3502028
#   IMF Primary Commodity Price System (PCPS):
#     https://www.imf.org/external/np/res/commod/index.asp
#
# =============================================================================

# -----------------------------------------------------------------------------
# 0. DEPENDENCY CHECK
# -----------------------------------------------------------------------------

if (!exists("combined")) {
  stop(paste(
    "Object 'combined' not found.",
    "Source 00_data_load_clean.R before running this script."
  ))
}

if (!exists("nb_model")) {
  message(paste(
    "Note: 'nb_model' not found in environment.",
    "Source 02_predictive.R for the baseline NB model to compare against.",
    "Continuing — extended model will still be fit."
  ))
}

# -----------------------------------------------------------------------------
# 1. PACKAGES
# -----------------------------------------------------------------------------

econ_packages <- c("wbstats", "countrycode", "dplyr", "tidyr", "ggplot2",
                   "MASS", "pscl", "car", "broom", "knitr", "kableExtra",
                   "scales", "readxl", "httr")

for (pkg in econ_packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg)
  }
}

suppressPackageStartupMessages({
  library(wbstats)
  library(countrycode)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(MASS)
  library(pscl)
  library(car)
  library(broom)
  library(knitr)
  library(kableExtra)
  library(scales)
})

if (!dir.exists("figures")) dir.create("figures")

# =============================================================================
# SECTION 1 — COUNTRY CROSSWALK
# =============================================================================
# SCAD country names do not always match World Bank names. The crosswalk
# below covers all non-obvious mismatches. The ~40 straightforward cases
# (where SCAD name == WB name) are handled by countrycode() for ISO codes.

wb_crosswalk <- c(
  "Democratic Republic of Congo"     = "Congo, Dem. Rep.",
  "Democratic Republic of the Congo" = "Congo, Dem. Rep.",
  "Republic of Congo"                = "Congo, Rep.",
  "Cote d'Ivoire"                    = "Cote d'Ivoire",
  "Ivory Coast"                      = "Cote d'Ivoire",
  "Egypt"                            = "Egypt, Arab Rep.",
  "Gambia"                           = "Gambia, The",
  "Swaziland"                        = "Eswatini",
  "Tanzania"                         = "Tanzania",
  "Sudan"                            = "Sudan",
  "South Sudan"                      = "South Sudan",
  "Trinidad and Tobago"              = "Trinidad and Tobago",
  "Dominican Republic"               = "Dominican Republic",
  "Yemen"                            = "Yemen, Rep."
)

combined <- combined %>%
  mutate(
    countryname_wb = if_else(
      countryname %in% names(wb_crosswalk),
      wb_crosswalk[countryname],
      countryname
    )
  )

combined <- combined %>%
  mutate(
    iso3c = countrycode(
      sourcevar   = countryname_wb,
      origin      = "country.name",
      destination = "iso3c",
      warn        = FALSE
    )
  )

scad_wb_crosswalk <- combined %>%
  distinct(countryname, countryname_wb, iso3c) %>%
  arrange(countryname)

# Exported scalars for Rmd inline reference
n_countries_crosswalk  <- nrow(scad_wb_crosswalk)
n_countries_matched    <- sum(!is.na(scad_wb_crosswalk$iso3c))
n_countries_unmatched  <- sum(is.na(scad_wb_crosswalk$iso3c))

if (n_countries_unmatched > 0) {
  message("Country crosswalk: ", n_countries_unmatched,
          " SCAD name(s) did not match an ISO3 code — check scad_wb_crosswalk.")
}

write.csv(scad_wb_crosswalk, "scad_wb_crosswalk.csv", row.names = FALSE)

# =============================================================================
# SECTION 2 — WORLD BANK GDP DATA
# =============================================================================

wb_indicators <- c(
  gdp_pc_const = "NY.GDP.PCAP.KD",
  gdp_growth   = "NY.GDP.MKTP.KD.ZG"
)

wb_raw <- tryCatch({
  wbstats::wb_data(
    indicator  = wb_indicators,
    country    = "all",
    start_date = 1990,
    end_date   = 2017
  )
}, error = function(e) {
  message("World Bank API call failed: ", conditionMessage(e),
          "\nCheck internet connection. Download WDI CSV manually from ",
          "https://databank.worldbank.org/source/world-development-indicators ",
          "if needed and load as a long-format data frame.")
  NULL
})

if (is.null(wb_raw)) {
  combined$gdp_pc_const <- NA_real_
  combined$gdp_growth   <- NA_real_
} else {
  wb_clean <- wb_raw %>%
    select(country, date, gdp_pc_const, gdp_growth) %>%
    rename(countryname_wb = country, styr = date) %>%
    filter(!is.na(styr)) %>%
    mutate(styr = as.integer(styr))

  combined <- combined %>%
    left_join(wb_clean, by = c("countryname_wb", "styr"))
}

# =============================================================================
# SECTION 3 — SHADOW ECONOMY DATA (Medina & Schneider 2019)
# =============================================================================
# Shadow economy = informal economic activity as % of official GDP.
#
# ACQUISITION STRATEGY (three routes, tried in order):
#   Route A: Direct download from Schneider's JKU academic page (Excel file)
#   Route B: Local IMF WP/18/17 appendix CSV (save as shadow_economy_imf2018.csv)
#   Route C: Regional average fallback — Africa ~39%, Latin America ~38%.
#            Flagged via shadow_proxy_used = TRUE; no country-level variation.

shadow_df <- NULL

if (requireNamespace("httr", quietly = TRUE) &&
    requireNamespace("readxl", quietly = TRUE)) {
  tryCatch({
    jku_excel_url <- paste0(
      "https://www.jku.at/fileadmin/gruppen/108/",
      "Research/Shadow_Economy/",
      "ShadowEconomies_Database.xlsx"
    )
    tmp <- tempfile(fileext = ".xlsx")
    resp <- httr::GET(jku_excel_url, httr::timeout(30),
                      httr::write_disk(tmp, overwrite = TRUE))
    if (httr::status_code(resp) == 200) {
      shadow_raw <- readxl::read_excel(tmp)
      shadow_df <- shadow_raw %>%
        rename_with(tolower) %>%
        select(matches("country|nation"), matches("year|date"), matches("shadow|informal")) %>%
        setNames(c("countryname_shadow", "styr", "shadow_pct")) %>%
        mutate(
          styr       = as.integer(styr),
          shadow_pct = as.numeric(shadow_pct)
        ) %>%
        filter(!is.na(styr), styr >= 1990, styr <= 2017)
    }
  }, error = function(e) {
    message("Shadow economy Route A failed: ", conditionMessage(e))
  })
}

if (is.null(shadow_df)) {
  imf_local_path <- "shadow_economy_imf2018.csv"
  if (file.exists(imf_local_path)) {
    shadow_df <- read.csv(imf_local_path, stringsAsFactors = FALSE) %>%
      mutate(
        styr       = as.integer(styr),
        shadow_pct = as.numeric(shadow_pct)
      ) %>%
      filter(!is.na(styr), styr >= 1990, styr <= 2017)
    # Note: IMF WP/18/17 covers 1991–2015 only; 2016–2017 will be NA.
  }
}

if (is.null(shadow_df)) {
  # Route C: regional average proxy (no country-level variation — document in A3)
  shadow_df <- combined %>%
    distinct(countryname_wb, styr, region) %>%
    mutate(
      shadow_pct = case_when(
        region == "Africa"        ~ 39.0,
        region == "Latin America" ~ 38.0,
        TRUE                      ~ 38.5
      ),
      countryname_shadow = countryname_wb
    ) %>%
    select(countryname_shadow, styr, shadow_pct)

  shadow_proxy_used <- TRUE
} else {
  shadow_proxy_used <- FALSE
}

if (!is.null(shadow_df) && !shadow_proxy_used) {
  shadow_df <- shadow_df %>%
    mutate(
      countryname_wb = case_when(
        grepl("Congo.*Dem|Dem.*Congo", countryname_shadow, ignore.case = TRUE) ~ "Congo, Dem. Rep.",
        grepl("Congo.*Rep|Republic.*Congo", countryname_shadow, ignore.case = TRUE) ~ "Congo, Rep.",
        grepl("Ivory|Cote d", countryname_shadow, ignore.case = TRUE) ~ "Cote d'Ivoire",
        grepl("Egypt", countryname_shadow, ignore.case = TRUE) ~ "Egypt, Arab Rep.",
        grepl("Gambia", countryname_shadow, ignore.case = TRUE) ~ "Gambia, The",
        grepl("Swaziland|Eswatini", countryname_shadow, ignore.case = TRUE) ~ "Eswatini",
        TRUE ~ countryname_shadow
      )
    )

  combined <- combined %>%
    left_join(
      shadow_df %>% select(countryname_wb, styr, shadow_pct),
      by = c("countryname_wb", "styr")
    )
} else if (shadow_proxy_used) {
  combined <- combined %>%
    left_join(
      shadow_df %>% select(countryname_shadow, styr, shadow_pct) %>%
        rename(countryname_wb = countryname_shadow),
      by = c("countryname_wb", "styr")
    )
}

combined <- combined %>%
  mutate(
    gdp_pc_adj     = gdp_pc_const * (1 + shadow_pct / 100),
    log_gdp_pc_adj = log(gdp_pc_adj)
  )

# =============================================================================
# SECTION 4 — IMF COMMODITY PRICE SHOCKS
# =============================================================================
# Each SCAD country is mapped to its primary commodity export.
# Commodity price shock = year-over-year % change in primary commodity price.
# Exogenous predictor: global commodity prices are not driven by any one country.
#
# ACQUISITION:
#   Download IMF commodity price CSV from:
#   https://www.imf.org/external/np/res/commod/index.asp
#   Save as "imf_commodity_prices.csv" in the working directory.
#   Fallback: World Bank commodity series via wbstats.

country_commodity_map <- list(
  "Nigeria"        = "oil",
  "Algeria"        = "oil",
  "Angola"         = "oil",
  "Libya"          = "oil",
  "Gabon"          = "oil",
  "Sudan"          = "oil",
  "South Sudan"    = "oil",
  "Chad"           = "oil",
  "Congo, Rep."    = "oil",
  "Cote d'Ivoire"  = "cocoa",
  "Cameroon"       = "cocoa",
  "Mali"           = "gold",
  "Burkina Faso"   = "gold",
  "Ghana"          = c("cocoa", "gold"),
  "Ethiopia"       = "coffee",
  "Uganda"         = "coffee",
  "Kenya"          = "coffee",
  "Rwanda"         = "coffee",
  "Tanzania"       = c("coffee", "gold"),
  "Congo, Dem. Rep." = "copper",
  "Zambia"           = "copper",
  "Niger"          = "uranium",
  "Egypt, Arab Rep." = "cotton",
  "Malawi"         = "tobacco",
  "Zimbabwe"       = "tobacco",
  "Guinea"         = "aluminum",
  "Cuba"              = "sugar",
  "Jamaica"           = "sugar",
  "Dominican Republic" = "sugar",
  "Haiti"             = "sugar",
  "Guatemala"      = "coffee",
  "Honduras"       = "coffee",
  "Nicaragua"      = "coffee",
  "El Salvador"    = "coffee"
)

imf_commodity_path <- "imf_commodity_prices.csv"
commodity_prices   <- NULL

if (file.exists(imf_commodity_path)) {
  raw_commodity <- read.csv(imf_commodity_path, stringsAsFactors = FALSE,
                             check.names = FALSE)
  commodity_prices <- raw_commodity %>%
    rename_with(tolower) %>%
    mutate(styr = as.integer(substr(as.character(date), 1, 4))) %>%
    filter(!is.na(styr), styr >= 1989, styr <= 2017) %>%
    group_by(styr) %>%
    summarise(across(where(is.numeric), mean, na.rm = TRUE), .groups = "drop")
}

if (is.null(commodity_prices)) {
  wb_commodity_raw <- tryCatch({
    wbstats::wb_data(
      indicator = c(
        oil_price    = "PNGASUSSDBBL",
        cocoa_price  = "PCOCOA",
        coffee_price = "PCOFFOTM",
        gold_price   = "PGOLD",
        copper_price = "PCOPP",
        agri_index   = "PAGRALL"
      ),
      country    = "WLD",
      start_date = 1989,
      end_date   = 2017
    )
  }, error = function(e) {
    message("WB commodity API failed: ", conditionMessage(e))
    NULL
  })

  if (!is.null(wb_commodity_raw)) {
    commodity_prices <- wb_commodity_raw %>%
      select(date, oil_price, cocoa_price, coffee_price,
             gold_price, copper_price, agri_index) %>%
      rename(styr = date) %>%
      mutate(styr = as.integer(styr))
  }
}

if (!is.null(commodity_prices)) {
  commodity_shocks <- commodity_prices %>%
    arrange(styr) %>%
    mutate(across(
      -styr,
      ~ (. - lag(.)) / lag(.) * 100,
      .names = "{.col}_shock"
    )) %>%
    select(styr, ends_with("_shock"))

  build_country_shocks <- function(country_name, commodities) {
    col_candidates <- paste0(commodities, "_price_shock")
    col_present    <- intersect(col_candidates, names(commodity_shocks))
    if (length(col_present) == 0) {
      col_present <- intersect("agri_index_shock", names(commodity_shocks))
    }
    commodity_shocks %>%
      select(styr, all_of(col_present)) %>%
      mutate(
        commodity_price_shock = rowMeans(select(., all_of(col_present)), na.rm = TRUE),
        countryname_wb        = country_name
      ) %>%
      select(countryname_wb, styr, commodity_price_shock)
  }

  country_shock_list <- mapply(
    build_country_shocks,
    country_name = names(country_commodity_map),
    commodities  = country_commodity_map,
    SIMPLIFY     = FALSE
  )
  country_shock_df <- bind_rows(country_shock_list)

  mapped_countries   <- names(country_commodity_map)
  all_wb_countries   <- unique(combined$countryname_wb)
  unmapped_countries <- setdiff(all_wb_countries, mapped_countries)

  if (length(unmapped_countries) > 0 && "agri_index_shock" %in% names(commodity_shocks)) {
    agri_fallback <- commodity_shocks %>%
      select(styr, agri_index_shock) %>%
      rename(commodity_price_shock = agri_index_shock) %>%
      crossing(countryname_wb = unmapped_countries)
    country_shock_df <- bind_rows(country_shock_df, agri_fallback)
  }

  combined <- combined %>%
    left_join(country_shock_df, by = c("countryname_wb", "styr"))
} else {
  combined$commodity_price_shock <- NA_real_
}

# =============================================================================
# SECTION 5 — LAG ECONOMIC PREDICTORS BY ONE YEAR
# =============================================================================
# GDP growth and commodity price shocks are lagged t-1 → t.
# GDP per capita level is a structural variable — no lag applied.

combined <- combined %>%
  group_by(countryname) %>%
  arrange(styr) %>%
  mutate(
    gdp_growth_lag1      = lag(gdp_growth,            1),
    commodity_shock_lag1 = lag(commodity_price_shock,  1)
  ) %>%
  ungroup()

scad_econ <- combined

# =============================================================================
# SECTION 6 — BUILD ECONOMIC MODEL DATASET
# =============================================================================

base_model_vars <- c(
  "ndeath", "high_sev",
  "etype", "escalation", "esc_high",
  "log_duration", "log_npart",
  "cgovtarget", "rgovtarget", "repress",
  "region", "styr_scaled",
  "cy_events_log", "countryname", "start_year"
)

if ("spill_hav" %in% names(scad_econ)) {
  spill_col_econ <- "spill_hav"
} else if ("spillover_context" %in% names(scad_econ)) {
  spill_col_econ <- "spillover_context"
} else {
  scad_econ$spillover_context <- 0L
  spill_col_econ <- "spillover_context"
}

econ_model_vars <- c(
  "log_gdp_pc_adj",
  "gdp_growth_lag1",
  "shadow_pct",
  "commodity_shock_lag1"
)

if (!"cy_events_log" %in% names(scad_econ)) {
  scad_econ <- scad_econ %>%
    group_by(countryname, start_year) %>%
    mutate(country_yr_events = n()) %>%
    ungroup() %>%
    mutate(cy_events_log = log1p(country_yr_events))
}

all_model_vars <- c(base_model_vars, spill_col_econ,
                    intersect(econ_model_vars, names(scad_econ)))

model_data_econ <- scad_econ %>%
  select(all_of(intersect(all_model_vars, names(scad_econ)))) %>%
  filter(
    !is.na(ndeath),
    !is.na(escalation),
    !is.na(log_duration),
    !is.na(log_npart),
    !is.na(log_gdp_pc_adj),
    !is.na(gdp_growth_lag1),
    !is.na(commodity_shock_lag1)
  )

model_data_econ$high_sev_f <- factor(
  model_data_econ$high_sev,
  levels = c(0, 1),
  labels = c("LowSev", "HighSev")
)

# Exported scalars for Rmd inline reference
n_model_econ    <- nrow(model_data_econ)
baseline_n_econ <- if (exists("model_df2")) nrow(model_df2) else
                   if (exists("model_df")) nrow(na.omit(model_df)) else NA
pct_retained_econ <- if (!is.na(baseline_n_econ) && baseline_n_econ > 0)
                       round(n_model_econ / baseline_n_econ * 100, 1) else NA
pct_lost_econ     <- if (!is.na(pct_retained_econ)) 100 - pct_retained_econ else NA

if (!is.na(pct_lost_econ) && pct_lost_econ > 20) {
  warning("Economic model: > 20% of baseline rows lost to economic data missingness (",
          pct_lost_econ, "%). Note as a limitation in the A3.")
}

# -- Table E1: Missingness by country ------------------------------------------

table_e1 <- scad_econ %>%
  group_by(countryname_wb, region) %>%
  summarise(
    n_events       = n(),
    miss_gdp_pc    = sum(is.na(gdp_pc_const)),
    miss_gdp_growth = sum(is.na(gdp_growth)),
    miss_shadow    = sum(is.na(shadow_pct)),
    miss_commod    = sum(is.na(commodity_price_shock)),
    .groups = "drop"
  ) %>%
  mutate(pct_miss_gdp = round(miss_gdp_pc / n_events * 100, 0)) %>%
  filter(pct_miss_gdp > 0 | miss_commod > 0) %>%
  arrange(desc(pct_miss_gdp))

kable(table_e1,
      col.names = c("Country (WB Name)", "Region", "Events",
                    "Miss GDP/Cap", "Miss GDP Growth",
                    "Miss Shadow", "Miss Commodity", "% Miss GDP"),
      caption   = "Table E1: Economic Data Missingness by Country")

# =============================================================================
# SECTION 7 — EXTENDED NEGATIVE BINOMIAL GLM
# =============================================================================

econ_predictors <- c("log_gdp_pc_adj", "gdp_growth_lag1", "commodity_shock_lag1")

if (!exists("shadow_proxy_used") || !shadow_proxy_used) {
  if ("shadow_pct" %in% names(model_data_econ) &&
      var(model_data_econ$shadow_pct, na.rm = TRUE) > 0) {
    econ_predictors <- c(econ_predictors, "shadow_pct")
  }
}

econ_terms_available <- intersect(econ_predictors, names(model_data_econ))

econ_formula_str <- paste(
  "ndeath ~",
  "factor(etype) + escalation + esc_high +",
  "log_duration + log_npart +",
  "cgovtarget + rgovtarget + repress +",
  "region + styr_scaled + cy_events_log +",
  spill_col_econ, "+",
  paste(econ_terms_available, collapse = " + ")
)

nb_formula_econ <- as.formula(econ_formula_str)

set.seed(2025)
train_idx_econ <- caret::createDataPartition(
  model_data_econ$high_sev_f, p = 0.75, list = FALSE
)
train_econ <- model_data_econ[ train_idx_econ, ]
test_econ  <- model_data_econ[-train_idx_econ, ]

nb_model_econ <- tryCatch(
  MASS::glm.nb(nb_formula_econ, data = train_econ),
  error = function(e) {
    warning("Extended NB GLM failed to converge: ", conditionMessage(e),
            "\nTry removing commodity_shock_lag1 first, then gdp_growth_lag1.")
    NULL
  }
)

if (!is.null(nb_model_econ)) {

  null_econ   <- MASS::glm.nb(ndeath ~ 1, data = train_econ)
  pseudo_econ <- round(1 - as.numeric(logLik(nb_model_econ)) /
                             as.numeric(logLik(null_econ)), 3)
  aic_econ    <- round(AIC(nb_model_econ), 1)

  econ_pred   <- predict(nb_model_econ, newdata = test_econ, type = "response")
  econ_corr   <- round(cor(log1p(test_econ$ndeath), log1p(econ_pred),
                            use = "complete.obs"), 3)
  econ_rmse   <- round(sqrt(mean((test_econ$ndeath - econ_pred)^2, na.rm = TRUE)), 2)

  assign("nb_model_econ", nb_model_econ, envir = .GlobalEnv)
  assign("pseudo_econ",   pseudo_econ,   envir = .GlobalEnv)
  assign("aic_econ",      aic_econ,      envir = .GlobalEnv)
  assign("econ_corr",     econ_corr,     envir = .GlobalEnv)
  assign("econ_rmse",     econ_rmse,     envir = .GlobalEnv)

  # -- Table E2: Coefficient comparison ----------------------------------------

  econ_tidy <- broom::tidy(nb_model_econ, conf.int = TRUE, exponentiate = FALSE) %>%
    filter(term != "(Intercept)") %>%
    mutate(
      term_clean = gsub("factor\\(etype\\)", "etype=", term),
      term_clean = gsub("regionLatin America", "Region: LatAm", term_clean),
      term_clean = gsub("_", " ", term_clean),
      coef_ci    = sprintf("%.3f [%.3f, %.3f]", estimate, conf.low, conf.high),
      sig_flag   = case_when(
        p.value < 0.001 ~ "***",
        p.value < 0.01  ~ "**",
        p.value < 0.05  ~ "*",
        p.value < 0.10  ~ ".",
        TRUE            ~ ""
      ),
      econ_var   = term %in% econ_terms_available
    ) %>%
    mutate(
      term_clean = if_else(econ_var, paste0("[ECON] ", term_clean), term_clean),
      p_rounded  = round(p.value, 4)
    ) %>%
    select(term_clean, coef_ci, p_rounded, sig_flag)

  kable(econ_tidy,
        col.names = c("Predictor", "Coef [95% CI]", "p-value", "Sig"),
        caption   = paste("Table E2: Extended NB GLM Coefficients.",
                          "[ECON] = new economic predictor.",
                          "Outcome: ndeath (negative binomial)."))

  # -- Table E3: GVIF -----------------------------------------------------------

  gvif_econ <- tryCatch({
    vif_raw <- car::vif(nb_model_econ, type = "predictor")
    if (is.matrix(vif_raw)) {
      data.frame(
        Predictor = rownames(vif_raw),
        GVIF      = round(vif_raw[, "GVIF"], 3),
        Df        = vif_raw[, "Df"],
        GVIF_adj  = round(vif_raw[, "GVIF^(1/(2*Df))"], 3)
      )
    } else {
      data.frame(
        Predictor = names(vif_raw),
        GVIF      = round(vif_raw, 3),
        Df        = 1,
        GVIF_adj  = round(sqrt(vif_raw), 3)
      )
    }
  }, error = function(e) {
    warning("GVIF computation failed: ", conditionMessage(e))
    NULL
  })

  if (!is.null(gvif_econ)) {
    gvif_econ <- gvif_econ %>%
      mutate(
        Flag = case_when(
          GVIF_adj > 3.16 ~ "HIGH (sqrt-GVIF > 3.16; investigate)",
          GVIF_adj > 2.24 ~ "MODERATE (sqrt-GVIF > 2.24; note)",
          TRUE            ~ "OK"
        )
      )

    kable(gvif_econ,
          col.names = c("Predictor", "GVIF", "Df", "sqrt-GVIF", "Flag"),
          caption   = paste("Table E3: GVIF Table — Extended NB GLM.",
                            "sqrt-GVIF thresholds: > 2.24 = moderate,",
                            "> 3.16 = high collinearity concern."))

    high_gvif <- gvif_econ %>% filter(grepl("HIGH", Flag)) %>% pull(Predictor)
    if (length(high_gvif) > 0) {
      warning("HIGH GVIF predictors: ", paste(high_gvif, collapse = ", "),
              ". Consider dropping or re-specifying. log_gdp_pc_adj and ",
              "gdp_growth_lag1 may be collinear.")
    }
  }

  # -- Table E4: Model performance comparison -----------------------------------

  table_e4 <- data.frame(
    Model = c("Baseline NB GLM", "Extended NB GLM (+ economic vars)"),
    N_train = c(
      if (exists("train_df")) nrow(train_df) else "N/A",
      nrow(train_econ)
    ),
    Pseudo_R2 = c(
      if (exists("nb_pseudo_r2")) nb_pseudo_r2 else "N/A",
      pseudo_econ
    ),
    AIC = c(
      if (exists("nb_aic")) nb_aic else "N/A",
      aic_econ
    ),
    Pearson_r_log = c(
      if (exists("nb_corr")) nb_corr else "N/A",
      econ_corr
    ),
    RMSE = c(
      if (exists("nb_rmse")) nb_rmse else "N/A",
      econ_rmse
    ),
    Meets_R2_Bound = c(
      if (exists("nb_pseudo_r2") && is.numeric(nb_pseudo_r2))
        ifelse(nb_pseudo_r2 >= 0.10, "YES", "NO") else "N/A",
      ifelse(pseudo_econ >= 0.10, "YES", "NO")
    )
  )

  kable(table_e4,
        col.names = c("Model", "Train N", "McFadden R²", "AIC",
                      "Pearson r (log)", "RMSE", "R² ≥ 0.10"),
        caption   = paste("Table E4: NB GLM Performance — Baseline vs. Extended.",
                          "Extended model uses economic model dataset (smaller N)."))
}

# =============================================================================
# SECTION 8 — ECONOMIC PREDICTOR COEFFICIENT PLOT (Figure E1)
# =============================================================================

if (!is.null(nb_model_econ)) {

  econ_coef_df <- broom::tidy(nb_model_econ, conf.int = TRUE,
                               exponentiate = FALSE) %>%
    filter(term %in% econ_terms_available) %>%
    mutate(
      term_label = case_when(
        term == "log_gdp_pc_adj"       ~ "Log GDP/capita\n(shadow-adjusted)",
        term == "gdp_growth_lag1"      ~ "GDP Growth Rate\n(%, lagged 1yr)",
        term == "commodity_shock_lag1" ~ "Commodity Price Shock\n(%, lagged 1yr)",
        term == "shadow_pct"           ~ "Shadow Economy\n(% of GDP)",
        TRUE                           ~ gsub("_", " ", term)
      ),
      significant = p.value < 0.05,
      direction   = ifelse(estimate > 0, "Positive", "Negative")
    )

  # Exported scalars for Rmd interpretation prose
  n_econ_predictors <- nrow(econ_coef_df)
  n_sig_econ        <- sum(econ_coef_df$significant, na.rm = TRUE)

  if (nrow(econ_coef_df) > 0) {
    fig_e1 <- ggplot(econ_coef_df,
                     aes(x = reorder(term_label, estimate),
                         y = estimate,
                         color = significant)) +
      geom_point(size = 3.5) +
      geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
                    width = 0.15, linewidth = 0.9) +
      geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
      scale_color_manual(
        values = c("TRUE" = "#C0392B", "FALSE" = "gray60"),
        labels = c("TRUE" = "p < 0.05", "FALSE" = "p ≥ 0.05"),
        name   = "Significance"
      ) +
      coord_flip() +
      labs(
        title    = "Figure E1: Economic Predictor Coefficients — Extended NB GLM",
        subtitle = paste0("Outcome: ndeath (negative binomial) | N = ",
                          nrow(train_econ), " training events"),
        x        = NULL,
        y        = "Coefficient (log-count scale)"
      ) +
      theme_minimal(base_size = 12) +
      theme(
        plot.title       = element_text(face = "bold"),
        panel.grid.minor = element_blank(),
        legend.position  = "bottom"
      )

    print(fig_e1)
    ggsave("figures/figE1_econ_coef_plot.png", fig_e1,
           width = 7, height = 5, dpi = 150)
  }
}

# =============================================================================
# SECTION 9 — UPDATED WHAT-IF COUNTERFACTUAL (Figure E2)
# =============================================================================
# Replicates the prescriptive counterfactual from 04_prescriptive.R
# (capping escalation at Level 3) using the extended NB GLM.
#
# CAVEATS (for Rmd prose, not cat()):
#   (1) In-model estimate only — not a causal guarantee.
#   (2) Capping escalation assumes intervention is feasible at Level 3.
#   (3) Economic model uses a subset of events; pct reduction applies to that subset.

if (!is.null(nb_model_econ) && nrow(test_econ) > 0) {

  test_econ_pred <- test_econ %>%
    mutate(pred_observed = predict(nb_model_econ, newdata = ., type = "response"))

  econ_obs_mean <- round(mean(test_econ_pred$pred_observed, na.rm = TRUE), 2)

  test_econ_cf <- test_econ %>%
    mutate(escalation = pmin(escalation, 3), esc_high = 0L)
  cf_preds_econ <- predict(nb_model_econ, newdata = test_econ_cf, type = "response")
  econ_cf_mean  <- round(mean(cf_preds_econ, na.rm = TRUE), 2)

  pct_reduction_econ <- round((econ_obs_mean - econ_cf_mean) / econ_obs_mean * 100, 1)

  # Baseline model comparison
  pct_reduction_base_econ <- NA_real_
  if (exists("nb_model") && exists("test_df")) {
    baseline_cf_e2 <- test_df %>%
      mutate(escalation = pmin(escalation, 3), esc_high = 0L)
    baseline_obs_mean_e2 <- mean(predict(nb_model, newdata = test_df,
                                          type = "response"), na.rm = TRUE)
    baseline_cf_mean_e2  <- mean(predict(nb_model, newdata = baseline_cf_e2,
                                          type = "response"), na.rm = TRUE)
    pct_reduction_base_econ <- round(
      (baseline_obs_mean_e2 - baseline_cf_mean_e2) / baseline_obs_mean_e2 * 100, 1
    )
  }

  # Exported scalars for Rmd inline reference
  assign("econ_obs_mean",            econ_obs_mean,            envir = .GlobalEnv)
  assign("econ_cf_mean",             econ_cf_mean,             envir = .GlobalEnv)
  assign("pct_reduction_econ",       pct_reduction_econ,       envir = .GlobalEnv)
  assign("pct_reduction_base_econ",  pct_reduction_base_econ,  envir = .GlobalEnv)

  cf_compare_e2 <- data.frame(
    Scenario = c("Observed\n(All Escalation Levels)",
                 "Counterfactual\n(Escalation Capped at Level 3)"),
    Mean_Deaths = c(econ_obs_mean, econ_cf_mean)
  )

  fig_e2 <- ggplot(cf_compare_e2,
                   aes(x = Scenario, y = Mean_Deaths, fill = Scenario)) +
    geom_col(width = 0.5, alpha = 0.85) +
    geom_text(aes(label = sprintf("%.2f deaths", Mean_Deaths)),
              vjust = -0.4, size = 4.5, fontface = "bold") +
    annotate("text", x = 1.5, y = max(cf_compare_e2$Mean_Deaths) * 0.85,
             label = sprintf("−%.1f%%\nreduction", pct_reduction_econ),
             size = 5, color = "#C0392B", fontface = "bold") +
    scale_fill_manual(values = c("#2C3E50", "#C0392B"), guide = "none") +
    scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
    labs(
      title    = "Figure E2: Counterfactual — Capping Escalation at Level 3",
      subtitle = paste0("Extended NB GLM with economic predictors | ",
                        "Test set N = ", nrow(test_econ)),
      x        = NULL,
      y        = "Mean Predicted Deaths per Event"
    ) +
    theme_minimal(base_size = 12) +
    theme(
      plot.title         = element_text(face = "bold"),
      panel.grid.minor   = element_blank(),
      panel.grid.major.x = element_blank()
    )

  print(fig_e2)
  ggsave("figures/figE2_whatif_counterfactual_econ.png", fig_e2,
         width = 7, height = 5, dpi = 150)
}

# =============================================================================
# SECTION 10 — UPDATED RELATIONAL CONNECTIONS (Additions to Table 7)
# =============================================================================

table_e5_additions <- data.frame(
  Utility = c(
    "log_gdp_pc_adj (shadow-adjusted GDP/capita)",
    "gdp_growth_lag1 (GDP growth rate, lagged 1yr)",
    "commodity_shock_lag1 (commodity price shock, lagged 1yr)"
  ),
  Analytical_Objective = c(
    "AO3: Economic impact on conflict severity (Extended NB GLM)",
    "AO3: Economic impact on conflict severity (Extended NB GLM)",
    "AO3: Economic impact on conflict severity (Extended NB GLM)"
  ),
  Managerial_Objective = c(
    "MO3: Predict conflict severity via economic indicators",
    "MO3: Predict conflict severity via economic indicators",
    "MO3: Predict conflict severity via economic indicators"
  ),
  Evidence = c(
    if (!is.null(nb_model_econ))
      paste0("Coef = ", round(coef(nb_model_econ)["log_gdp_pc_adj"], 3),
             "; p = ", round(summary(nb_model_econ)$coef["log_gdp_pc_adj", 4], 4))
    else "Extended model not available",
    if (!is.null(nb_model_econ) && "gdp_growth_lag1" %in% names(coef(nb_model_econ)))
      paste0("Coef = ", round(coef(nb_model_econ)["gdp_growth_lag1"], 3),
             "; p = ", round(summary(nb_model_econ)$coef["gdp_growth_lag1", 4], 4))
    else "Extended model not available",
    if (!is.null(nb_model_econ) && "commodity_shock_lag1" %in% names(coef(nb_model_econ)))
      paste0("Coef = ", round(coef(nb_model_econ)["commodity_shock_lag1"], 3),
             "; p = ", round(summary(nb_model_econ)$coef["commodity_shock_lag1", 4], 4))
    else "Extended model not available"
  ),
  Current_Slack = c(
    "GDP data available for ~80-90% of events; fragile states are the gap",
    "Lagged growth accounts for endogeneity; pre-conflict economic conditions",
    "Exogenous price shock variable; avoids reverse causality from conflict"
  )
)

kable(table_e5_additions,
      col.names = c("Utility (Predictor)", "Analytical Objective",
                    "Managerial Objective", "Evidence", "Current Slack"),
      caption   = "Table E5: New Connections — Economic Utilities to MO3")