# =============================================================================
# 00_data_load_clean.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
# Data: SCAD v3.3 (Social Conflict Analysis Database) — public domain
# =============================================================================
#
# PURPOSE:
#   This is the foundational script. Every other script (01–04) depends on
#   the objects created here. Run this first, or source() it from the master
#   Rmd setup chunk before any other source() calls.
#
# OUTPUT OBJECTS (available in global environment after sourcing):
#   africa_raw       — raw Africa CSV as loaded
#   latam_raw        — raw Latin America CSV as loaded
#   africa           — cleaned Africa data with all derived variables
#   latam            — cleaned Latin America data with all derived variables
#   combined         — row-bound africa + latam with region indicator
#   model_df         — combined, complete-cases only, analysis-ready
#
# ENCODING NOTE:
#   Both CSVs require latin1 encoding. UTF-8 fails due to special characters
#   in countryname and source note fields. Do not change this.
#
# SENTINEL VALUES:
#   SCAD uses -99 and -88 as coded missing values (not actual negatives).
#   These are replaced with NA throughout before any analysis.
#   Key fields affected: ndeath, npart, duration, escalation.
#
# DERIVED VARIABLES CREATED:
#   high_sev         — binary flag: ndeath > 10; NA when ndeath is NA
#   log_ndeath       — log1p(ndeath), stabilizes right-skewed distribution
#   log_duration     — log1p(duration), removes right-skew for regression
#   log_npart        — log1p(npart), same rationale
#   etype_f          — etype as labeled factor (11 levels per SCAD codebook)
#   esc_high         — binary: escalation >= 4 (threshold effect term)
#   styr_scaled      — year centered on 1990 (removes multicollinearity with
#                      intercept; coefficient = deaths per additional year)
#   lat_bin / lon_bin — 2.5-degree grid cell bins for coarse spillover proxy
#   region           — "Africa" or "Latin America" character column
#
# =============================================================================

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

required_packages <- c(
  "dplyr", "tidyr", "lubridate", "ggplot2", "scales",
  "MASS", "pscl", "randomForest", "caret", "pROC",
  "geosphere", "corrplot", "car", "broom",
  "knitr", "kableExtra", "gridExtra"
)

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

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(lubridate)
  library(ggplot2)
  library(scales)
})

# -----------------------------------------------------------------------------
# 2. FILE PATHS
# -----------------------------------------------------------------------------
# *** CHANGE THESE TO YOUR LOCAL PATHS BEFORE RUNNING ***
# All team members: update these two lines only. Everything else is portable.

AFRICA_PATH <- "SCAD2018Africa_Final.csv"
LATAM_PATH  <- "SCAD2018LatinAmerica_Final.csv"

# Assumes the project directory where the CSVs are co-located; point these
# at your own local copies if running from elsewhere.

# -----------------------------------------------------------------------------
# 3. LOAD RAW DATA
# -----------------------------------------------------------------------------

africa_raw <- tryCatch(
  read.csv(AFRICA_PATH, encoding = "latin1", stringsAsFactors = FALSE),
  error = function(e) stop(paste(
    "Could not load Africa CSV.\n",
    "Check AFRICA_PATH at the top of 00_data_load_clean.R.\n",
    "Error:", conditionMessage(e)
  ))
)

latam_raw <- tryCatch(
  read.csv(LATAM_PATH, encoding = "latin1", stringsAsFactors = FALSE),
  error = function(e) stop(paste(
    "Could not load Latin America CSV.\n",
    "Check LATAM_PATH at the top of 00_data_load_clean.R.\n",
    "Error:", conditionMessage(e)
  ))
)

# -----------------------------------------------------------------------------
# 4. CLEANING FUNCTION
# -----------------------------------------------------------------------------

clean_region <- function(df, region_label) {

  # -- 4a. Tag region ----------------------------------------------------------
  df$region <- region_label

  # -- 4b. Replace sentinel missing values with NA -----------------------------
  # SCAD uses -99 and -88 as coded missing. These are NOT real values.
  sentinel_fields <- c("ndeath", "npart", "duration", "escalation",
                       "latitude", "longitude")
  for (col in sentinel_fields) {
    if (col %in% names(df)) {
      df[[col]] <- ifelse(df[[col]] %in% c(-99, -88), NA, df[[col]])
    }
  }

  # Also replace any remaining negatives in ndeath with NA
  if ("ndeath" %in% names(df)) {
    df$ndeath <- ifelse(!is.na(df$ndeath) & df$ndeath < 0, NA, df$ndeath)
  }

  # -- 4c. Parse dates ---------------------------------------------------------
  if ("startdate" %in% names(df)) {
    df$startdate <- as.Date(df$startdate, format = "%Y-%m-%d")
  }
  if ("enddate" %in% names(df)) {
    df$enddate <- as.Date(df$enddate, format = "%Y-%m-%d")
  }

  if ("startdate" %in% names(df)) {
    df$start_year  <- as.integer(format(df$startdate, "%Y"))
    df$start_month <- as.integer(format(df$startdate, "%m"))
  } else if ("styr" %in% names(df)) {
    df$start_year  <- df$styr
    df$start_month <- df$stmo
  }

  # -- 4d. Filter to 1990+ for comparability ----------------------------------
  if ("start_year" %in% names(df)) {
    df <- df[!is.na(df$start_year) & df$start_year >= 1990, ]
  }

  # -- 4e. Derived variables ---------------------------------------------------

  # Binary high-severity outcome — NA when ndeath is NA (not silently 0)
  df$high_sev <- ifelse(
    is.na(df$ndeath),
    NA_integer_,
    ifelse(df$ndeath > 10, 1L, 0L)
  )

  df$log_ndeath   <- log1p(pmax(df$ndeath,   0, na.rm = FALSE))
  df$log_duration <- log1p(pmax(df$duration, 0, na.rm = FALSE))
  df$log_npart    <- log1p(pmax(df$npart,    0, na.rm = FALSE))

  # Escalation threshold binary: Level 4+ marks qualitative jump in severity
  df$esc_high <- ifelse(!is.na(df$escalation) & df$escalation >= 4, 1L, 0L)

  # Year scaled to 1990 baseline
  if ("start_year" %in% names(df)) {
    df$styr_scaled <- df$start_year - 1990
  }

  # Event type as labeled factor — 11 levels per SCAD v3.3 codebook
  # etype=11 "Mob Violence" added in SCAD v3.3; treat as NA if absent in data
  etype_labels <- c(
    "1"  = "Government Violence",
    "2"  = "Demonstration",
    "3"  = "Organized Violence",
    "4"  = "Spontaneous Violence",
    "5"  = "Strike",
    "6"  = "Lockout",
    "7"  = "Seizure",
    "8"  = "Mutiny",
    "9"  = "Riot",
    "10" = "Disruptive State Act",
    "11" = "Mob Violence"
  )
  if ("etype" %in% names(df)) {
    df$etype_f <- factor(df$etype,
                         levels = as.integer(names(etype_labels)),
                         labels = unname(etype_labels))
  }

  # 2.5-degree grid cell bins (coarse proxy — Haversine flag computed in 03)
  if (all(c("latitude", "longitude") %in% names(df))) {
    df$lat_bin <- floor(df$latitude  / 2.5) * 2.5
    df$lon_bin <- floor(df$longitude / 2.5) * 2.5
  }

  # -- 4f. Recode government targeting to clean binary -------------------------
  if ("cgovtarget" %in% names(df)) {
    df$cgovtarget <- ifelse(is.na(df$cgovtarget), 0L, as.integer(df$cgovtarget == 1))
  }
  if ("rgovtarget" %in% names(df)) {
    df$rgovtarget <- ifelse(is.na(df$rgovtarget), 0L, as.integer(df$rgovtarget == 1))
  }
  if ("repress" %in% names(df)) {
    df$repress <- ifelse(is.na(df$repress), 0L, as.integer(df$repress == 1))
  }

  return(df)
}

# -----------------------------------------------------------------------------
# 5. APPLY CLEANING
# -----------------------------------------------------------------------------

africa <- clean_region(africa_raw, "Africa")
latam  <- clean_region(latam_raw,  "Latin America")

# -----------------------------------------------------------------------------
# 6. COMBINE REGIONS
# -----------------------------------------------------------------------------

combined <- bind_rows(africa, latam)
combined$region <- factor(combined$region, levels = c("Africa", "Latin America"))

# -----------------------------------------------------------------------------
# 7. MODEL-READY DATASET
# -----------------------------------------------------------------------------

model_vars <- c(
  "ndeath", "log_ndeath", "high_sev",
  "etype", "escalation", "esc_high",
  "log_duration", "log_npart",
  "cgovtarget", "rgovtarget", "repress",
  "region", "styr_scaled",
  "latitude", "longitude",
  "lat_bin", "lon_bin",
  "start_year", "start_month",
  "countryname"
)

model_vars_present <- intersect(model_vars, names(combined))

model_df <- combined[, model_vars_present] %>%
  filter(
    !is.na(ndeath),
    !is.na(escalation),
    !is.na(log_duration),
    !is.na(log_npart),
    !is.na(latitude),
    !is.na(longitude)
  )

# -----------------------------------------------------------------------------
# 8. VALIDATION SCALARS (available for inline reference in Rmd)
# -----------------------------------------------------------------------------

n_africa_model    <- sum(model_df$region == "Africa")
n_latam_model     <- sum(model_df$region == "Latin America")
n_model_df        <- nrow(model_df)
n_dropped         <- nrow(combined) - nrow(model_df)
neg_deaths_remain <- sum(model_df$ndeath < 0, na.rm = TRUE)
hs_rate_model     <- round(mean(model_df$high_sev, na.rm = TRUE) * 100, 1)