# =============================================================================
# 01_descriptive.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:
#   Descriptive stage of the IIA Descriptive → Predictive → Prescriptive
#   pipeline. Produces all summary tables, statistical tests, and figures
#   for the "Measure / Current State" section of the S-A3.
#
# DEPENDS ON:
#   00_data_load_clean.R must be sourced first.
#   Objects required: africa, latam, combined, model_df
#
# OUTPUTS:
#   Table 1  — Regional summary statistics
#   Table 2  — Statistical tests (t-tests, chi-square)
#   Figure 1 — Event type distribution by region
#   Figure 2 — Mean deaths by event type (log scale)
#   Figure 3 — Annual event counts and death totals 1990–2017
#   Figure 4 — Monthly seasonality by region
#   Figure 5 — Geographic dot map (4,000-event sample, opacity = deaths)
#   Figure 6 — Mean deaths by escalation level (log scale)
#
# =============================================================================

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

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(scales)
  library(knitr)
  library(kableExtra)
  library(gridExtra)
})

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

theme_conflict <- function() {
  theme_minimal(base_size = 12) +
    theme(
      plot.title    = element_text(face = "bold", size = 13),
      plot.subtitle = element_text(color = "gray40", size = 10),
      legend.position = "bottom",
      panel.grid.minor = element_blank(),
      strip.text = element_text(face = "bold")
    )
}

region_colors <- c("Africa" = "#C0392B", "Latin America" = "#2980B9")

# =============================================================================
# TABLE 1: REGIONAL SUMMARY STATISTICS
# =============================================================================

regional_summary <- combined %>%
  group_by(region) %>%
  summarise(
    Events             = n(),
    Total_Deaths       = sum(ndeath, na.rm = TRUE),
    Mean_Deaths        = round(mean(ndeath, na.rm = TRUE), 2),
    Median_Deaths      = median(ndeath, na.rm = TRUE),
    Max_Deaths         = max(ndeath, na.rm = TRUE),
    Mean_Duration_days = round(mean(duration, na.rm = TRUE), 1),
    HighSev_N          = sum(high_sev, na.rm = TRUE),
    HighSev_Pct        = round(mean(high_sev, na.rm = TRUE) * 100, 1),
    .groups = "drop"
  )

combined_row <- combined %>%
  summarise(
    region             = factor("Combined"),
    Events             = n(),
    Total_Deaths       = sum(ndeath, na.rm = TRUE),
    Mean_Deaths        = round(mean(ndeath, na.rm = TRUE), 2),
    Median_Deaths      = median(ndeath, na.rm = TRUE),
    Max_Deaths         = max(ndeath, na.rm = TRUE),
    Mean_Duration_days = round(mean(duration, na.rm = TRUE), 1),
    HighSev_N          = sum(high_sev, na.rm = TRUE),
    HighSev_Pct        = round(mean(high_sev, na.rm = TRUE) * 100, 1)
  )

table1 <- bind_rows(regional_summary, combined_row)

kable(table1,
      col.names = c("Region", "Events", "Total Deaths", "Mean Deaths",
                    "Median Deaths", "Max Deaths", "Mean Duration (days)",
                    "High-Sev N", "High-Sev %"),
      caption = "Table 1: Regional Summary Statistics — SCAD v3.3 (1990–2017)")

# Scalars referenced in Rmd inline text
africa_mean_deaths <- regional_summary$Mean_Deaths[regional_summary$region == "Africa"]
latam_mean_deaths  <- regional_summary$Mean_Deaths[regional_summary$region == "Latin America"]
mortality_ratio    <- round(africa_mean_deaths / latam_mean_deaths, 1)

# =============================================================================
# TABLE 2: STATISTICAL TESTS
# =============================================================================

af_deaths  <- africa$ndeath[!is.na(africa$ndeath) & africa$ndeath >= 0]
la_deaths  <- latam$ndeath[!is.na(latam$ndeath)   & latam$ndeath  >= 0]
af_dur     <- africa$duration[!is.na(africa$duration) & africa$duration >= 0]
la_dur     <- latam$duration[!is.na(latam$duration)  & latam$duration  >= 0]

ttest_deaths   <- t.test(af_deaths, la_deaths)
ttest_duration <- t.test(af_dur, la_dur)

hs_table   <- table(combined$region, combined$high_sev)
chisq_hs   <- chisq.test(hs_table)

etype_table  <- table(combined$region, combined$etype)
chisq_etype  <- chisq.test(etype_table)

table2 <- data.frame(
  Test = c(
    "t-test: Mean deaths (Africa vs. LatAm)",
    "t-test: Mean duration (Africa vs. LatAm)",
    "Chi-square: High-severity rate by region",
    "Chi-square: Event type distribution by region"
  ),
  Statistic = c(
    round(ttest_deaths$statistic, 3),
    round(ttest_duration$statistic, 3),
    round(chisq_hs$statistic, 3),
    round(chisq_etype$statistic, 3)
  ),
  p_value = c(
    format(ttest_deaths$p.value,   scientific = TRUE, digits = 3),
    format(ttest_duration$p.value, scientific = TRUE, digits = 3),
    format(chisq_hs$p.value,       scientific = TRUE, digits = 3),
    format(chisq_etype$p.value,    scientific = TRUE, digits = 3)
  ),
  Significant = c(
    ttest_deaths$p.value   < 0.05,
    ttest_duration$p.value < 0.05,
    chisq_hs$p.value       < 0.05,
    chisq_etype$p.value    < 0.05
  )
)

kable(table2, caption = "Table 2: Regional Differences — Statistical Tests")

# =============================================================================
# FIGURE 1: EVENT TYPE DISTRIBUTION BY REGION
# =============================================================================

etype_dist <- combined %>%
  filter(!is.na(etype_f)) %>%
  count(region, etype_f) %>%
  group_by(region) %>%
  mutate(pct = n / sum(n) * 100) %>%
  ungroup()

fig1 <- ggplot(etype_dist, aes(x = reorder(etype_f, pct), y = pct, fill = region)) +
  geom_col(position = "dodge", alpha = 0.85, width = 0.7) +
  coord_flip() +
  scale_fill_manual(values = region_colors) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  labs(
    title    = "Figure 1: Event Type Distribution by Region",
    subtitle = "SCAD v3.3 (1990–2017) — percentage of total events per region",
    x        = NULL,
    y        = "% of Regional Events",
    fill     = NULL
  ) +
  theme_conflict()

print(fig1)
ggsave("figures/fig1_etype_distribution.png", fig1, width = 8, height = 5, dpi = 150)

# =============================================================================
# FIGURE 2: MEAN DEATHS BY EVENT TYPE (LOG SCALE)
# =============================================================================
# Key finding: etype 1 (Government Violence) and etype 3 (Organized Violence)
# are the deadliest despite NOT being the most frequent event types.

etype_deaths <- combined %>%
  filter(!is.na(etype_f), !is.na(ndeath), ndeath >= 0) %>%
  group_by(region, etype_f) %>%
  summarise(
    mean_deaths = mean(ndeath, na.rm = TRUE),
    n_events    = n(),
    .groups = "drop"
  ) %>%
  filter(n_events >= 10)

fig2 <- ggplot(etype_deaths,
               aes(x = reorder(etype_f, mean_deaths), y = mean_deaths, fill = region)) +
  geom_col(position = "dodge", alpha = 0.85, width = 0.7) +
  coord_flip() +
  scale_y_log10(labels = comma) +
  scale_fill_manual(values = region_colors) +
  labs(
    title    = "Figure 2: Mean Estimated Deaths by Event Type",
    subtitle = "Log scale — Government & Organized Violence are deadliest despite lower frequency",
    x        = NULL,
    y        = "Mean Deaths per Event (log scale)",
    fill     = NULL
  ) +
  theme_conflict()

print(fig2)
ggsave("figures/fig2_deaths_by_etype.png", fig2, width = 8, height = 5, dpi = 150)

# =============================================================================
# FIGURE 3: ANNUAL EVENT COUNTS AND DEATH TOTALS (1990–2017)
# =============================================================================

annual_trends <- combined %>%
  filter(!is.na(start_year), start_year >= 1990, start_year <= 2017) %>%
  group_by(region, start_year) %>%
  summarise(
    n_events     = n(),
    total_deaths = sum(ndeath, na.rm = TRUE),
    .groups = "drop"
  )

fig3a <- ggplot(annual_trends, aes(x = start_year, y = n_events,
                                    color = region, group = region)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
  scale_color_manual(values = region_colors) +
  scale_y_continuous(labels = comma) +
  scale_x_continuous(breaks = seq(1990, 2017, by = 5)) +
  labs(
    title    = "Figure 3a: Annual Conflict Events (1990–2017)",
    subtitle = "Africa shows accelerating trend post-2010 (Arab Spring / Sahel instability)",
    x        = "Year",
    y        = "Number of Events",
    color    = NULL
  ) +
  theme_conflict()

fig3b <- ggplot(annual_trends, aes(x = start_year, y = total_deaths,
                                    color = region, group = region)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.5) +
  scale_color_manual(values = region_colors) +
  scale_y_continuous(labels = comma) +
  scale_x_continuous(breaks = seq(1990, 2017, by = 5)) +
  labs(
    title    = "Figure 3b: Annual Estimated Deaths (1990–2017)",
    subtitle = "Death totals are episodic; driven by a small number of mass-casualty events",
    x        = "Year",
    y        = "Total Estimated Deaths",
    color    = NULL
  ) +
  theme_conflict()

grid.arrange(fig3a, fig3b, ncol = 1)
ggsave("figures/fig3a_annual_events.png", fig3a, width = 8, height = 4, dpi = 150)
ggsave("figures/fig3b_annual_deaths.png", fig3b, width = 8, height = 4, dpi = 150)

# =============================================================================
# FIGURE 4: MONTHLY SEASONALITY
# =============================================================================
# Africa shows a secondary July peak not present in Latin America.
# Consistent with dry-season agricultural conflict literature (sub-Saharan).

month_labels <- c("Jan","Feb","Mar","Apr","May","Jun",
                  "Jul","Aug","Sep","Oct","Nov","Dec")

monthly <- combined %>%
  filter(!is.na(start_month)) %>%
  count(region, start_month) %>%
  group_by(region) %>%
  mutate(
    pct        = n / sum(n) * 100,
    month_abbr = factor(month_labels[start_month], levels = month_labels)
  ) %>%
  ungroup()

fig4 <- ggplot(monthly, aes(x = month_abbr, y = pct,
                             color = region, group = region)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2.5) +
  scale_color_manual(values = region_colors) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  labs(
    title    = "Figure 4: Seasonal Conflict Frequency by Month",
    subtitle = "Both regions peak Jan–Feb; Africa shows secondary July peak (dry-season pattern)",
    x        = "Month",
    y        = "% of Annual Events",
    color    = NULL
  ) +
  theme_conflict()

print(fig4)
ggsave("figures/fig4_seasonality.png", fig4, width = 8, height = 4, dpi = 150)

# =============================================================================
# FIGURE 5: GEOGRAPHIC DOT MAP
# =============================================================================
# 4,000-event sample; opacity scaled to ndeath so high-severity events
# are visually prominent. Geocoding artifact (capital-city clustering)
# addressed in 03_spillover.R choropleth (Fig 5A).

set.seed(2025)
map_sample <- combined %>%
  filter(!is.na(latitude), !is.na(longitude),
         !is.na(ndeath), ndeath >= 0) %>%
  sample_n(min(4000, n()))

death_cap  <- quantile(map_sample$ndeath, 0.99, na.rm = TRUE)
map_sample <- map_sample %>%
  mutate(death_alpha = pmin(ndeath, death_cap) / death_cap)

fig5 <- ggplot(map_sample, aes(x = longitude, y = latitude,
                                color = region, alpha = death_alpha)) +
  geom_point(size = 0.8) +
  scale_color_manual(values = region_colors) +
  scale_alpha_continuous(range = c(0.05, 0.85), guide = "none") +
  coord_quickmap() +
  labs(
    title    = "Figure 5: Conflict Event Locations (4,000-event sample)",
    subtitle = "Point opacity proportional to estimated deaths — darker = more lethal",
    x        = "Longitude",
    y        = "Latitude",
    color    = NULL,
    caption  = "Note: Geocoding artifacts cluster some events at capital cities.\nSee choropleth (Fig 5A) in supplementary analysis for corrected view."
  ) +
  theme_conflict() +
  theme(plot.caption = element_text(color = "gray50", size = 8))

print(fig5)
ggsave("figures/fig5_geographic_map.png", fig5, width = 9, height = 5, dpi = 150)

# =============================================================================
# FIGURE 6: MEAN DEATHS BY ESCALATION LEVEL
# =============================================================================
# Mean deaths inflect sharply at escalation Level 4 — the "intervention window."
# This is the empirical basis for the esc_high binary term in predictive models.

esc_summary <- combined %>%
  filter(!is.na(escalation), !is.na(ndeath), ndeath >= 0) %>%
  group_by(region, escalation) %>%
  summarise(
    mean_deaths = mean(ndeath, na.rm = TRUE),
    n_events    = n(),
    .groups = "drop"
  ) %>%
  filter(n_events >= 5)

fig6 <- ggplot(esc_summary,
               aes(x = factor(escalation), y = mean_deaths,
                   color = region, group = region)) +
  geom_line(linewidth = 1.1) +
  geom_point(aes(size = n_events), alpha = 0.8) +
  scale_color_manual(values = region_colors) +
  scale_y_log10(labels = comma) +
  scale_size_continuous(range = c(2, 7), guide = "none") +
  geom_vline(xintercept = which(levels(factor(esc_summary$escalation)) == "4"),
             linetype = "dashed", color = "gray40", linewidth = 0.8) +
  annotate("text", x = which(levels(factor(esc_summary$escalation)) == "4") + 0.15,
           y = max(esc_summary$mean_deaths, na.rm = TRUE) * 0.6,
           label = "Level 4\ninflection", size = 3, color = "gray30", hjust = 0) +
  labs(
    title    = "Figure 6: Mean Estimated Deaths by Escalation Level",
    subtitle = "Log scale — sharp increase at Level 4 marks the intervention threshold",
    x        = "Escalation Level",
    y        = "Mean Deaths per Event (log scale)",
    color    = NULL
  ) +
  theme_conflict()

print(fig6)
ggsave("figures/fig6_escalation_deaths.png", fig6, width = 8, height = 5, dpi = 150)

# =============================================================================
# SUPPLEMENTARY: GOVERNMENT TARGETING
# =============================================================================

gov_target_summary <- combined %>%
  group_by(region) %>%
  summarise(
    n_events       = n(),
    cgovtarget_n   = sum(cgovtarget, na.rm = TRUE),
    cgovtarget_pct = round(mean(cgovtarget, na.rm = TRUE) * 100, 1),
    rgovtarget_n   = sum(rgovtarget, na.rm = TRUE),
    rgovtarget_pct = round(mean(rgovtarget, na.rm = TRUE) * 100, 1),
    .groups = "drop"
  )

kable(gov_target_summary,
      col.names = c("Region", "Events", "Central Gov Target N",
                    "Central Gov Target %", "Regional Gov Target N",
                    "Regional Gov Target %"),
      caption   = "Supplementary: Government-Targeted Events by Region")

# =============================================================================
# SUPPLEMENTARY: TOP COUNTRIES
# =============================================================================

top_countries <- combined %>%
  filter(!is.na(countryname)) %>%
  count(region, countryname, sort = TRUE) %>%
  group_by(region) %>%
  slice_head(n = 5) %>%
  mutate(rank = row_number()) %>%
  ungroup()

kable(top_countries,
      col.names = c("Region", "Country", "Events", "Rank"),
      caption   = "Supplementary: Top 5 Countries by Conflict Event Count")