---
title: "Economic Hardship in Arizona"
subtitle: "A Community Analytics Dashboard — PAF 516 Final Project"
format:
  dashboard:
    orientation: columns
    embed-resources: true
execute:
  echo: false
  warning: false
  message: false
---

```{r}
#| label: setup
#| include: false
options(tigris_use_cache = TRUE)
options(tigris_progress_bar = FALSE)
set.seed(42)
```

```{r}
#| label: packages
#| include: false
if (!file.exists("~/.paf516_ready")) {
  if (!requireNamespace("renv", quietly = TRUE)) install.packages("renv")
  tmp <- tempfile(fileext = ".lock")
  download.file(
    url      = "https://raw.githubusercontent.com/AntJam-Howell/paf-516-labs/main/renv.lock",
    destfile = tmp, mode = "wb"
  )
  renv::restore(lockfile = tmp, prompt = FALSE)
  file.create("~/.paf516_ready")
}
library(tidyverse); library(tidycensus); library(sf); library(spdep)
library(viridis); library(patchwork); library(tigris); library(DT)
library(mapgl); library(glue); library(plotly); library(scales)
```

```{r}
#| label: county-config
# ── STUDENT: change these 3 lines + update the YAML subtitle above ──────
# Examples:
#   Los Angeles:      TARGET_STATE="CA", TARGET_COUNTY="Los Angeles",  TARGET_LABEL="Los Angeles County, CA"
#   Cook (Chicago):   TARGET_STATE="IL", TARGET_COUNTY="Cook",         TARGET_LABEL="Cook County, IL"
#   Miami-Dade:       TARGET_STATE="FL", TARGET_COUNTY="Miami-Dade",   TARGET_LABEL="Miami-Dade County, FL"
#   King (Seattle):   TARGET_STATE="WA", TARGET_COUNTY="King",         TARGET_LABEL="King County, WA"
#   Harris (Houston): TARGET_STATE="TX", TARGET_COUNTY="Harris",       TARGET_LABEL="Harris County, TX"
TARGET_STATE  <- "AZ"
TARGET_COUNTY <- "Maricopa"
TARGET_LABEL  <- "Maricopa County, AZ"
# ─────────────────────────────────────────────────────────────────────────
```

```{r}
#| label: index-config
# ╔══════════════════════════════════════════════════════════════════════════╗
# ║              STUDENT CONFIGURATION — READ CAREFULLY                     ║
# ╚══════════════════════════════════════════════════════════════════════════╝
#
# The BASELINE Economic Hardship Index (EHI) uses THREE components:
#   1. Poverty rate          (higher = more hardship)
#   2. Unemployment rate     (higher = more hardship)
#   3. Median household income, inverted (lower income = more hardship)
#
# YOUR TASK: Expand the index by uncommenting ONE or TWO options below.
# After editing, re-render the dashboard. Then compare how the spatial
# patterns, cluster maps, and policy findings shift with your new index.
#
# ── AVAILABLE EXTRA COMPONENTS ──────────────────────────────────────────────
#
#  Option A  Renter Burden
#            % households that rent (housing cost pressure)
#            ACS: B25003_003 (renter occ.) / B25003_001 (total occ.)
#
#  Option B  Low Educational Attainment
#            % adults with no high school diploma
#            ACS: B06009_002 (< HS) / B06009_001 (pop. 25+)
#
#  Option C  Food Insecurity (SNAP receipt)
#            % households receiving SNAP/food stamps
#            ACS: B22003_002 (received SNAP) / B22003_001 (total HH)
#
#  Option D  Transportation Disadvantage
#            % households with no vehicle available
#            ACS: B08201_002 (0 vehicles) / B08201_001 (total HH)
#
# ── SET YOUR CHOICE HERE ──────────────────────────────────────────────────
# NOTE: After changing STUDENT_COMPONENTS, delete the _cache/ folder next to
# this file (or use Session > Restart R) before re-rendering, so that data
# pulls refresh with your new variables.
# Uncomment the component(s) you want to add.

STUDENT_COMPONENTS <- c(
#   "renter_rate",       # Option A — Renter Burden
#   "no_hs_rate",        # Option B — Low Educational Attainment
#   "snap_rate",         # Option C — Food Insecurity (SNAP)
#  "no_vehicle_rate"    # Option D — Transportation Disadvantage  ← ACTIVE
)

# ── INTERNALS — do not edit below this line ──────────────────────────────
EXTRA_VAR_MENU <- list(
  renter_rate = list(
    vars      = c(renter_num = "B25003_003", renter_den = "B25003_001"),
    rate_fn   = function(df) df$renter_num / df$renter_den,
    direction = +1,
    label     = "Renter Burden"
  ),
  no_hs_rate = list(
    vars      = c(no_hs_num = "B06009_002", no_hs_den = "B06009_001"),
    rate_fn   = function(df) df$no_hs_num / df$no_hs_den,
    direction = +1,
    label     = "Low Ed. Attainment"
  ),
  snap_rate = list(
    vars      = c(snap_num = "B22003_002", snap_den = "B22003_001"),
    rate_fn   = function(df) df$snap_num / df$snap_den,
    direction = +1,
    label     = "Food Insecurity (SNAP)"
  ),
  no_vehicle_rate = list(
    vars      = c(no_veh_num = "B08201_002", no_veh_den = "B08201_001"),
    rate_fn   = function(df) df$no_veh_num / df$no_veh_den,
    direction = +1,
    label     = "Transp. Disadvantage"
  )
)

extra_vars <- if (length(STUDENT_COMPONENTS) > 0)
  unlist(lapply(STUDENT_COMPONENTS, function(c) EXTRA_VAR_MENU[[c]]$vars)) else c()

hardship_vars <- c(
  pov_below50   = "C17002_002",
  pov_50to99    = "C17002_003",
  pov_den       = "C17002_001",
  unemp_num     = "B23025_005",
  unemp_den     = "B23025_002",
  median_income = "B19013_001"
)
all_hardship_vars <- c(hardship_vars, extra_vars)

INDEX_N     <- 3 + length(STUDENT_COMPONENTS)
INDEX_LABEL <- paste(
  c("Poverty", "Unemployment", "Income (inv.)",
    unlist(lapply(STUDENT_COMPONENTS, function(c) EXTRA_VAR_MENU[[c]]$label))),
  collapse = " + "
)

build_hardship_index <- function(raw_df, extra_components = STUDENT_COMPONENTS) {
  df <- raw_df %>%
    select(-moe) %>%
    pivot_wider(names_from = variable, values_from = estimate) %>%
    mutate(
      poverty_rate = (pov_below50 + pov_50to99) / pov_den,
      unemp_rate   = unemp_num / unemp_den
    ) %>%
    filter(!is.na(poverty_rate), !is.na(unemp_rate), !is.na(median_income))

  for (comp in extra_components)
    df[[comp]] <- EXTRA_VAR_MENU[[comp]]$rate_fn(df)

  (function(df) {
    geom <- if (inherits(df, "sf")) st_geometry(df) else NULL
    tbl  <- st_drop_geometry(df)
    z_mat <- cbind(
      as.numeric(scale(tbl$poverty_rate)),
      as.numeric(scale(tbl$unemp_rate)),
      -1 * as.numeric(scale(tbl$median_income))
    )
    for (comp in extra_components) {
      d <- EXTRA_VAR_MENU[[comp]]$direction
      z_mat <- cbind(z_mat, d * as.numeric(scale(tbl[[comp]])))
    }
    tbl$hardship_index <- rowMeans(z_mat, na.rm = TRUE)
    if (!is.null(geom)) st_sf(tbl, geometry = geom) else tbl
  })(df)
}
```

```{r}
#| label: pull-national
#| cache: true
county_raw <- get_acs(
  geography    = "county",
  variables    = all_hardship_vars,
  year         = 2023,
  survey       = "acs5",
  geometry     = TRUE,
  progress_bar = FALSE
)
county_2023 <- build_hardship_index(county_raw) %>%
  filter(!str_sub(GEOID, 1, 2) %in% c("72","78","66","60","69"))
```

```{r}
#| label: pull-az-tracts
#| cache: true
#| cache.extra: !expr paste(sort(STUDENT_COMPONENTS), collapse = ",")
# All Arizona census tracts — state-level interactive map + AZ LISA
az_tract_raw <- get_acs(
  geography    = "tract",
  variables    = all_hardship_vars,
  state        = TARGET_STATE,
  year         = 2023,
  survey       = "acs5",
  geometry     = TRUE,
  progress_bar = FALSE
)
az_tract_2023 <- build_hardship_index(az_tract_raw)

# Maricopa County boundary — yellow outline on AZ maps
maricopa_boundary <- az_tract_2023 %>%
  filter(str_detect(NAME, TARGET_COUNTY)) %>%
  st_union() %>%
  st_sf() %>%
  st_transform(4326)
```

```{r}
#| label: lisa-az
#| cache: true
#| cache.extra: !expr paste(sort(STUDENT_COMPONENTS), collapse = ",")
# AZ-wide LISA — 2023 cross-sectional hardship clustering across all AZ tracts
az_valid <- az_tract_2023 %>%
  filter(!is.na(hardship_index))

nb_az  <- poly2nb(az_valid, queen = TRUE)
w_az   <- nb2listw(nb_az, style = "W", zero.policy = TRUE)
lisa_az_res <- localmoran(az_valid$hardship_index, w_az, zero.policy = TRUE)

z_az   <- scale(az_valid$hardship_index)[, 1]
wz_az  <- lag.listw(w_az, z_az, zero.policy = TRUE)
p_az   <- lisa_az_res[, "Pr(z != E(Ii))"]

az_valid <- az_valid %>%
  mutate(
    lisa_az = case_when(
      z_az > 0 & wz_az > 0 & p_az < 0.05 ~ "HH",
      z_az < 0 & wz_az < 0 & p_az < 0.05 ~ "LL",
      z_az > 0 & wz_az < 0 & p_az < 0.05 ~ "HL",
      z_az < 0 & wz_az > 0 & p_az < 0.05 ~ "LH",
      TRUE                                 ~ "NS"
    )
  )

moran_az   <- moran.test(az_valid$hardship_index, w_az, zero.policy = TRUE)
az_n_hh    <- sum(az_valid$lisa_az == "HH", na.rm = TRUE)
az_n_ll    <- sum(az_valid$lisa_az == "LL", na.rm = TRUE)
az_moran_i <- round(unname(moran_az$estimate[1]), 3)

# Prepare for mapgl — simplify + transform + ensure character column
az_lisa_map <- az_valid %>%
  mutate(
    lisa_az       = factor(lisa_az, levels = c("HH","LL","HL","LH","NS")),
    lisa_az_label = case_when(
      lisa_az == "HH" ~ "Hot Spot \u2014 High Hardship Cluster",
      lisa_az == "LL" ~ "Cold Spot \u2014 Low Hardship Cluster",
      lisa_az == "HL" ~ "High-Low Outlier",
      lisa_az == "LH" ~ "Low-High Outlier",
      TRUE            ~ "Not Significant"
    ),
    tooltip_text = paste0(
      str_extract(NAME, "Census Tract [\\d\\.]+"), "\n",
      lisa_az_label, "\n",
      "Economic Hardship Index: ", round(hardship_index, 3)
    )
  ) %>%
  st_simplify(preserveTopology = TRUE, dTolerance = 500) %>%
  st_transform(4326)
```

```{r}
#| label: pull-local-current
#| cache: true
#| cache.extra: !expr paste(sort(STUDENT_COMPONENTS), collapse = ",")
local_combined_raw <- get_acs(
  geography    = "tract",
  variables    = all_hardship_vars,
  state        = TARGET_STATE,
  county       = TARGET_COUNTY,
  year         = 2023,
  survey       = "acs5",
  geometry     = TRUE,
  progress_bar = FALSE
)
local_2023 <- build_hardship_index(local_combined_raw)
```

```{r}
#| label: pull-temporal
#| cache: true
#| cache.extra: !expr paste(sort(STUDENT_COMPONENTS), collapse = ",")
local_2013_raw <- get_acs(
  geography = "tract", variables = all_hardship_vars,
  state = TARGET_STATE, county = TARGET_COUNTY,
  year = 2013, survey = "acs5", geometry = TRUE, progress_bar = FALSE
)
local_2019_raw <- get_acs(
  geography = "tract", variables = all_hardship_vars,
  state = TARGET_STATE, county = TARGET_COUNTY,
  year = 2019, survey = "acs5", geometry = TRUE, progress_bar = FALSE
)
local_2016_raw <- get_acs(
  geography = "tract", variables = all_hardship_vars,
  state = TARGET_STATE, county = TARGET_COUNTY,
  year = 2016, survey = "acs5", geometry = FALSE, progress_bar = FALSE
)
```

```{r}
#| label: pooled-std
# Extended to include student-selected extra components in pooled standardization
widen_and_rate <- function(df) {
  out <- df %>%
    select(-moe) %>%
    pivot_wider(names_from = variable, values_from = estimate) %>%
    mutate(
      poverty_rate = (pov_below50 + pov_50to99) / pov_den,
      unemp_rate   = unemp_num / unemp_den
    ) %>%
    filter(!is.na(poverty_rate), !is.na(unemp_rate), !is.na(median_income))
  for (comp in STUDENT_COMPONENTS)
    out[[comp]] <- EXTRA_VAR_MENU[[comp]]$rate_fn(out)
  out
}

local_2013_wide <- widen_and_rate(local_2013_raw)
local_2019_wide <- widen_and_rate(local_2019_raw)

vars_pool <- c("poverty_rate", "unemp_rate", "median_income", STUDENT_COMPONENTS)
r2013 <- local_2013_wide %>% st_drop_geometry() %>% select(GEOID, all_of(vars_pool))
r2019 <- local_2019_wide %>% st_drop_geometry() %>% select(GEOID, all_of(vars_pool))

pooled_stats <- map_dfr(vars_pool, function(v) {
  vals <- c(r2013[[v]], r2019[[v]])
  tibble(variable = v, mean = mean(vals, na.rm = TRUE), sd = sd(vals, na.rm = TRUE))
})

apply_pooled_z <- function(df, stats) {
  for (v in c("poverty_rate", "unemp_rate", "median_income")) {
    df[[paste0("z_", v)]] <- (df[[v]] - stats$mean[stats$variable == v]) /
                               stats$sd[stats$variable == v]
  }
  df$z_median_income <- -1 * df$z_median_income
  for (comp in STUDENT_COMPONENTS) {
    d <- EXTRA_VAR_MENU[[comp]]$direction
    df[[paste0("z_", comp)]] <- d *
      (df[[comp]] - stats$mean[stats$variable == comp]) /
       stats$sd[stats$variable == comp]
  }
  extra_z <- if (length(STUDENT_COMPONENTS) > 0) paste0("z_", STUDENT_COMPONENTS) else character(0)
  z_cols <- c("z_poverty_rate", "z_unemp_rate", "z_median_income", extra_z)
  z_tbl  <- st_drop_geometry(df)[, z_cols, drop = FALSE]
  df$hardship_index <- rowMeans(as.matrix(z_tbl), na.rm = TRUE)
  df
}

r2013_z <- apply_pooled_z(r2013, pooled_stats)
r2019_z <- apply_pooled_z(r2019, pooled_stats)

change_df <- local_2019_wide %>%
  select(GEOID, NAME) %>%
  left_join(r2019_z %>% select(GEOID, hardship_2019 = hardship_index), by = "GEOID") %>%
  inner_join(r2013_z %>% select(GEOID, hardship_2013 = hardship_index), by = "GEOID") %>%
  mutate(hardship_change = hardship_2019 - hardship_2013)

lo      <- quantile(change_df$hardship_change, 0.02, na.rm = TRUE)
hi      <- quantile(change_df$hardship_change, 0.98, na.rm = TRUE)
max_abs <- max(abs(lo), abs(hi))
change_df <- change_df %>%
  mutate(change_plot = pmax(pmin(hardship_change, max_abs), -max_abs))
```

```{r}
#| label: lisa-analysis
# Maricopa-only LISA — used for trajectory analysis (2019 pooled z-scores)
nb   <- poly2nb(change_df, queen = TRUE)
w    <- nb2listw(nb, style = "W", zero.policy = TRUE)
lisa <- localmoran(change_df$hardship_2019, w, zero.policy = TRUE)
z19  <- scale(change_df$hardship_2019)[,1]
wz19 <- lag.listw(w, z19, zero.policy = TRUE)
p19  <- lisa[, "Pr(z != E(Ii))"]

change_df <- change_df %>%
  mutate(lisa_2019 = case_when(
    z19 > 0 & wz19 > 0 & p19 < 0.05 ~ "HH",
    z19 < 0 & wz19 < 0 & p19 < 0.05 ~ "LL",
    z19 > 0 & wz19 < 0 & p19 < 0.05 ~ "HL",
    z19 < 0 & wz19 > 0 & p19 < 0.05 ~ "LH",
    TRUE ~ "NS"
  ))
```

```{r}
#| label: trajectories
lisa_2013 <- localmoran(change_df$hardship_2013, w, zero.policy = TRUE)
z13  <- scale(change_df$hardship_2013)[,1]
wz13 <- lag.listw(w, z13, zero.policy = TRUE)
p13  <- lisa_2013[, "Pr(z != E(Ii))"]

local_2016_wide <- widen_and_rate(local_2016_raw %>% mutate(geometry = NULL))
r2016_z <- apply_pooled_z(
  local_2016_wide %>% select(GEOID, all_of(vars_pool)),
  pooled_stats
)

change_df <- change_df %>%
  mutate(
    lisa_2013 = case_when(
      z13 > 0 & wz13 > 0 & p13 < 0.05 ~ "HH",
      z13 < 0 & wz13 < 0 & p13 < 0.05 ~ "LL",
      z13 > 0 & wz13 < 0 & p13 < 0.05 ~ "HL",
      z13 < 0 & wz13 > 0 & p13 < 0.05 ~ "LH",
      TRUE ~ "NS"
    ),
    trajectory = case_when(
      lisa_2013 == "HH" & lisa_2019 == "HH" ~ "Persistent HH",
      lisa_2013 != "HH" & lisa_2019 == "HH" ~ "Emerging HH",
      lisa_2013 == "HH" & lisa_2019 != "HH" ~ "Dissolving HH",
      lisa_2013 == "LL" & lisa_2019 == "LL" ~ "Persistent LL",
      lisa_2013 != "LL" & lisa_2019 == "LL" ~ "Emerging LL",
      lisa_2013 == "HL" | lisa_2019 == "HL" ~ "HL Outlier",
      TRUE                                   ~ "Stable NS"
    )
  ) %>%
  left_join(r2016_z %>% select(GEOID, hardship_2016 = hardship_index), by = "GEOID") %>%
  mutate(trend_type = case_when(
    hardship_2013 > hardship_2016 & hardship_2016 > hardship_2019 ~ "Consistently improving",
    hardship_2013 < hardship_2016 & hardship_2016 < hardship_2019 ~ "Consistently worsening",
    hardship_2013 > hardship_2016 & hardship_2016 < hardship_2019 ~ "V-shaped (improved mid)",
    hardship_2013 < hardship_2016 & hardship_2016 > hardship_2019 ~ "Inverted-V (worsened mid)",
    TRUE ~ "Flat / non-monotone"
  ))
```

```{r}
#| label: summary-stats
moran_global <- moran.test(change_df$hardship_2019, w, zero.policy = TRUE)

# Maricopa trajectory stats (for Tab 3 + Policy)
n_hh           <- sum(change_df$lisa_2019 == "HH", na.rm = TRUE)
n_ll           <- sum(change_df$lisa_2019 == "LL", na.rm = TRUE)
n_tracts       <- nrow(change_df)
pct_improved   <- round(100 * mean(change_df$hardship_change < -0.05, na.rm=TRUE), 1)
pct_worsened   <- round(100 * mean(change_df$hardship_change >  0.05, na.rm=TRUE), 1)
pct_persist_hh <- round(100 * mean(change_df$trajectory == "Persistent HH", na.rm=TRUE), 1)
pct_emerg_hh   <- round(100 * mean(change_df$trajectory == "Emerging HH",   na.rm=TRUE), 1)
n_persist_hh   <- sum(change_df$trajectory == "Persistent HH", na.rm = TRUE)
n_emerg_hh     <- sum(change_df$trajectory == "Emerging HH",   na.rm = TRUE)

n_counties <- sum(!is.na(county_2023$hardship_index))

traj_colors <- c(
  "Persistent HH" = "#db2b27", "Emerging HH"   = "#fdbf11",
  "Dissolving HH" = "#fdd870", "Persistent LL" = "#1696d2",
  "Emerging LL"   = "#73bfe2", "HL Outlier"    = "#ec008b",
  "Stable NS"     = "#d2d2d2"
)
lisa_colors <- c(
  "HH" = "#db2b27", "LL" = "#1696d2",
  "HL" = "#fdbf11", "LH" = "#73bfe2", "NS" = "#d2d2d2"
)

theme_urbn_dash <- function(base_size = 11) {
  theme_void(base_size = base_size) +
  theme(
    plot.background   = element_rect(fill = "#ffffff", color = NA),
    panel.background  = element_rect(fill = "#ffffff", color = NA),
    plot.title        = element_text(color = "#333333", size = base_size + 1,
                                     face = "bold", margin = margin(b = 3)),
    plot.subtitle     = element_text(color = "#767676", size = base_size - 1,
                                     margin = margin(b = 6)),
    plot.caption      = element_text(color = "#767676", size = base_size - 2,
                                     hjust = 0, margin = margin(t = 4)),
    legend.background = element_rect(fill = "#ffffff", color = NA),
    legend.key        = element_rect(fill = "#ffffff", color = NA),
    legend.text       = element_text(color = "#333333", size = base_size - 1),
    legend.title      = element_text(color = "#767676", size = base_size - 1),
    plot.margin       = margin(8, 8, 8, 8)
  )
}

# National context rankings
us_ranks <- county_2023 %>%
  st_drop_geometry() %>%
  filter(!is.na(hardship_index)) %>%
  mutate(us_rank = rank(-hardship_index, ties.method = "min")) %>%
  select(GEOID, us_rank)

az_counties <- county_2023 %>%
  st_drop_geometry() %>%
  filter(str_detect(NAME, "Arizona")) %>%
  left_join(us_ranks, by = "GEOID") %>%
  arrange(desc(hardship_index)) %>%
  mutate(
    az_rank     = row_number(),
    n_az        = n(),
    county_name = str_remove(NAME, " County, Arizona"),
    County      = county_name,
    `EH Index Pos. = Worse-Off` = round(hardship_index, 3),
    `AZ Rank (1 = most hardship, out of 15)` = az_rank,
    `US Rank (1 = most hardship, out of ~3,140)` = us_rank
  ) %>%
  select(County,
         `EH Index Pos. = Worse-Off`,
         `AZ Rank (1 = most hardship, out of 15)`,
         `US Rank (1 = most hardship, out of ~3,140)`)

# Within-county spatial heterogeneity table (all 15 AZ counties)
az_county_hetero <- az_tract_2023 %>%
  st_drop_geometry() %>%
  filter(!is.na(hardship_index)) %>%
  mutate(
    county_name = str_extract(NAME, "(?<=, )[^,]+(?=, Arizona)"),
    county_name = str_remove(county_name, " County")
  ) %>%
  group_by(county_name) %>%
  summarise(
    n_tracts = n(),
    mean_hi  = round(mean(hardship_index,   na.rm = TRUE), 3),
    sd_hi    = round(sd(hardship_index,     na.rm = TRUE), 3),
    median_hi= round(median(hardship_index, na.rm = TRUE), 3),
    min_hi   = round(min(hardship_index,    na.rm = TRUE), 3),
    max_hi   = round(max(hardship_index,    na.rm = TRUE), 3),
    .groups  = "drop"
  ) %>%
  arrange(desc(sd_hi)) %>%
  rename(
    County          = county_name,
    `N Tracts`      = n_tracts,
    `Mean`          = mean_hi,
    `SD (↑ = more variation)` = sd_hi,
    `Median`        = median_hi,
    `Min`           = min_hi,
    `Max`           = max_hi
  )
```

# National Context

## Column {width=55%}

```{r}
#| label: map-national
#| title: "Economic Hardship Across U.S. Counties"
county_map <- county_2023 %>%
  filter(!is.na(hardship_index)) %>%
  st_simplify(preserveTopology = TRUE, dTolerance = 5000) %>%
  st_transform(4326)

az_outline <- states(cb = TRUE) %>%
  filter(STUSPS == TARGET_STATE) %>%
  st_transform(4326)

maplibre(
  style  = carto_style("dark-matter"),
  bounds = c(-125.0, 24.0, -66.0, 50.0)
) |>
  add_fill_layer(
    id = "hardship-fill",
    source = county_map,
    fill_color = interpolate(
      column = "hardship_index",
      values = c(-2, -1, 0, 1, 2),
      stops  = c("#1696d2", "#73bfe2", "#f5f5f5", "#ec6a5a", "#db2b27")
    ),
    fill_opacity = 0.8,
    tooltip = "NAME"
  ) |>
  add_line_layer(
    id = "county-borders",
    source = county_map,
    line_color = "#333333",
    line_width = 0.2
  ) |>
  add_line_layer(
    id = "az-highlight",
    source = az_outline,
    line_color = "#fdbf11",
    line_width = 2.5
  ) |>
  htmlwidgets::onRender("
    function(el, x) {
      var leg = document.createElement('div');
      leg.style.cssText = 'position:absolute;bottom:32px;left:10px;background:rgba(15,15,15,0.82);padding:7px 10px 6px;border-radius:4px;font-family:Lato,sans-serif;min-width:180px;';
      leg.innerHTML = '<div style=\"font-size:10px;color:#cccccc;font-weight:600;margin-bottom:5px;\">Economic Hardship Index</div>'
        + '<div style=\"display:flex;align-items:center;gap:6px;\">'
        +   '<span style=\"font-size:9px;color:#73bfe2;white-space:nowrap;\">Better-off</span>'
        +   '<div style=\"flex:1;height:10px;border-radius:2px;background:linear-gradient(to right,#1696d2,#73bfe2,#f5f5f5,#ec6a5a,#db2b27);\"></div>'
        +   '<span style=\"font-size:9px;color:#ec6a5a;white-space:nowrap;\">Worse-off</span>'
        + '</div>';
      el.appendChild(leg);
    }
  ")
```

## Column {width=45%}

### Row {height=35%}

```{r}
#| label: table-az-counties
#| title: "Economic Hardship (EH) Rankings — Arizona Counties"
n_us_counties <- nrow(us_ranks)

az_counties %>%
  DT::datatable(
    colnames = c(
      "County",
      "EH Index",
      "AZ Rank",
      "US Rank"
    ),
    options = list(
      pageLength = 8, dom = "tp",
      columnDefs = list(list(className = "dt-center", targets = 1:3)),
      headerCallback = DT::JS(
        "function(thead) {",
        "  $(thead).find('th').eq(2).html('AZ Rank<br><small style=\"font-weight:normal;color:#777;\">out of 15 Counties</small>');",
        paste0("  $(thead).find('th').eq(3).html('US Rank<br><small style=\"font-weight:normal;color:#777;\">out of ", n_us_counties, " Counties</small>');"),
        "}"
      )
    ),
    rownames = FALSE
  ) %>%
  DT::formatStyle(
    columns = "County",
    target  = "row",
    backgroundColor = DT::styleEqual(TARGET_COUNTY, "#fffde6")
  )
```

### Row {height=65%}

```{r}
#| label: card-dot-compare
#| title: "All Arizona Counties vs. U.S. Extremes"
top5_us <- county_2023 %>%
  st_drop_geometry() %>%
  left_join(us_ranks, by = "GEOID") %>%
  filter(!is.na(hardship_index)) %>%
  slice_min(us_rank, n = 5, with_ties = FALSE) %>%
  mutate(
    state_name   = str_extract(NAME, "(?<=, ).+"),
    state_abb    = state.abb[match(state_name, state.name)],
    county_short = str_remove(NAME, " County.*"),
    group = "Top 5 US (Most Economic Hardship)",
    label = paste0(county_short, ", ", coalesce(state_abb, state_name))
  )

bot5_us <- county_2023 %>%
  st_drop_geometry() %>%
  left_join(us_ranks, by = "GEOID") %>%
  filter(!is.na(hardship_index)) %>%
  slice_max(us_rank, n = 5, with_ties = FALSE) %>%
  mutate(
    state_name   = str_extract(NAME, "(?<=, ).+"),
    state_abb    = state.abb[match(state_name, state.name)],
    county_short = str_remove(NAME, " County.*"),
    group = "Bottom 5 US (Least Economic Hardship)",
    label = paste0(county_short, ", ", coalesce(state_abb, state_name))
  )

az_all <- county_2023 %>%
  st_drop_geometry() %>%
  left_join(us_ranks, by = "GEOID") %>%
  filter(str_detect(NAME, "Arizona"), !is.na(hardship_index)) %>%
  mutate(
    group = "Arizona Counties",
    label = paste0(str_remove(NAME, " County, Arizona"), ", AZ")
  )

compare_df <- bind_rows(top5_us, az_all, bot5_us) %>%
  mutate(
    color_group = case_when(
      group == "Arizona Counties"                    ~ "#fdbf11",
      group == "Top 5 US (Most Economic Hardship)"   ~ "#db2b27",
      TRUE                                           ~ "#1696d2"
    ),
    label = fct_reorder(label, hardship_index)
  )

az_mean   <- mean(compare_df$hardship_index[compare_df$group == "Arizona Counties"],
                  na.rm = TRUE)
mid_label <- levels(compare_df$label)[ceiling(nlevels(compare_df$label) / 2)]

p <- ggplot(compare_df,
       aes(x = hardship_index, y = label, color = color_group,
           text = paste0(label, "\nIndex: ", round(hardship_index, 3)))) +
  geom_point(size = 2.8) +
  scale_color_identity() +
  theme_urbn_dash(base_size = 10) +
  theme(
    axis.text.y        = element_text(size = 8, color = "#333333"),
    axis.text.x        = element_text(size = 9, color = "#767676"),
    panel.grid.major.x = element_line(color = "#eeeeee", linewidth = 0.4)
  ) +
  labs(x = "Economic Hardship Index", y = NULL)

ggplotly(p, tooltip = "text") %>%
  layout(
    title  = list(
      text = paste0("All Arizona Counties vs. U.S. Extremes",
                    "<br><sup><span style=\"color:#fdbf11\">&#9679;</span> AZ Counties &nbsp;",
                    "<span style=\"color:#db2b27\">&#9679;</span> 5 Highest US &nbsp;",
                    "<span style=\"color:#1696d2\">&#9679;</span> 5 Lowest US &nbsp;",
                    "<span style=\"color:#9d9d9d\">&#9474;</span> AZ Mean</sup>"),
      font = list(size = 12, color = "#333333")
    ),
    margin = list(t = 50, b = 30),
    xaxis  = list(title = list(text = "Economic Hardship Index",
                               font = list(size = 11, color = "#555555"))),
    yaxis  = list(automargin = TRUE),
    shapes = list(list(
      type = "line",
      x0 = az_mean, x1 = az_mean, y0 = 0, y1 = 1, yref = "paper",
      line = list(color = "#9d9d9d", width = 1.5, dash = "dot")
    ))
  ) %>%
  config(displayModeBar = FALSE)
```

# Arizona in Focus {orientation="rows"}

## Top Row {height=55%}

### EHI Map {width=33%}

```{r}
#| label: map-az-tracts
#| title: "Economic Hardship Index: All Arizona Census Tracts (2023)"
az_map <- az_tract_2023 %>%
  filter(!is.na(hardship_index)) %>%
  st_simplify(preserveTopology = TRUE, dTolerance = 500) %>%
  st_transform(4326)

maplibre(
  style  = carto_style("dark-matter"),
  bounds = c(-115.2, 31.0, -108.8, 37.5)
) |>
  add_fill_layer(
    id           = "az-tracts",
    source       = az_map,
    fill_color   = interpolate(
      column = "hardship_index",
      values = c(-2, -1, 0, 1, 2),
      stops  = c("#1696d2", "#73bfe2", "#f5f5f5", "#ec6a5a", "#db2b27")
    ),
    fill_opacity = 0.8,
    tooltip      = "NAME"
  ) |>
  add_line_layer(
    id           = "maricopa-border",
    source       = maricopa_boundary,
    line_color   = "#fdbf11",
    line_width   = 2.5
  ) |>
  htmlwidgets::onRender("
    function(el, x) {
      // Maricopa boundary badge
      var badge = document.createElement('div');
      badge.style.cssText = 'position:absolute;top:10px;left:10px;background:rgba(10,10,10,0.75);color:#fdbf11;font-size:11px;font-family:Lato,sans-serif;padding:5px 9px;border-radius:3px;border:1px solid #fdbf11;pointer-events:none;';
      badge.innerHTML = '&#9632;&nbsp;Maricopa County boundary';
      el.appendChild(badge);
      // Compact gradient legend
      var leg = document.createElement('div');
      leg.style.cssText = 'position:absolute;bottom:32px;left:10px;background:rgba(15,15,15,0.82);padding:7px 10px 6px;border-radius:4px;font-family:Lato,sans-serif;min-width:170px;';
      leg.innerHTML = '<div style=\"font-size:10px;color:#cccccc;font-weight:600;margin-bottom:5px;\">Economic Hardship Index</div>'
        + '<div style=\"display:flex;align-items:center;gap:6px;\">'
        +   '<span style=\"font-size:9px;color:#73bfe2;white-space:nowrap;\">Better-off</span>'
        +   '<div style=\"flex:1;height:10px;border-radius:2px;background:linear-gradient(to right,#1696d2,#73bfe2,#f5f5f5,#ec6a5a,#db2b27);\"></div>'
        +   '<span style=\"font-size:9px;color:#ec6a5a;white-space:nowrap;\">Worse-off</span>'
        + '</div>';
      el.appendChild(leg);
    }
  ")
```

### Decomposing the Economic Hardship Index {width=34%}

```{r}
#| label: scatter-pov-unemp
#| title: "Decomposing the Economic Hardship Index"
county_name_lookup_scatter <- county_2023 %>%
  st_drop_geometry() %>%
  filter(str_detect(NAME, "Arizona")) %>%
  mutate(
    county_geoid = GEOID,
    county_name  = str_remove(str_extract(NAME, "^[^,]+"), " County")
  ) %>%
  select(county_geoid, county_name)

scatter_df <- az_tract_2023 %>%
  st_drop_geometry() %>%
  filter(!is.na(poverty_rate), !is.na(unemp_rate), !is.na(hardship_index)) %>%
  mutate(county_geoid = substr(GEOID, 1, 5)) %>%
  left_join(county_name_lookup_scatter, by = "county_geoid") %>%
  mutate(
    is_maricopa = county_name == TARGET_COUNTY,
    tract_label = paste0(
      str_extract(NAME, "Census Tract [\\d\\.]+"), "\n",
      coalesce(county_name, ""), " County\n",
      "Poverty: ",    scales::percent(poverty_rate, accuracy = 0.1), "\n",
      "Unemployed: ", scales::percent(unemp_rate,   accuracy = 0.1), "\n",
      "EHI: ",        round(hardship_index, 3)
    )
  )

# Compute 2D kernel density per tract to modulate alpha:
# dense clusters stay vivid; sparse outliers (high pov + high unemp) fade
dens <- MASS::kde2d(scatter_df$poverty_rate, scatter_df$unemp_rate, n = 100,
                    lims = c(range(scatter_df$poverty_rate), range(scatter_df$unemp_rate)))
get_density <- function(x, y, dens) {
  xi <- findInterval(x, dens$x)
  yi <- findInterval(y, dens$y)
  xi <- pmax(1L, pmin(xi, length(dens$x) - 1L))
  yi <- pmax(1L, pmin(yi, length(dens$y) - 1L))
  dens$z[cbind(xi, yi)]
}
scatter_df <- scatter_df %>%
  mutate(
    pt_density = get_density(poverty_rate, unemp_rate, dens),
    pt_alpha   = scales::rescale(pt_density, to = c(0.72, 0.92))
  )

p_scatter <- ggplot(scatter_df,
    aes(x = poverty_rate, y = unemp_rate, color = hardship_index,
        text = tract_label, size = is_maricopa, alpha = pt_alpha)) +
  geom_point() +
  scale_color_gradient2(
    low = "#1696d2", mid = "#2d2d2d", high = "#db2b27",
    midpoint = 0, name = "EHI", limits = c(-3, 3),
    oob = scales::squish       # squish OOB values to nearest limit color, not NA/gray
  ) +
  scale_size_manual(values = c("FALSE" = 1.2, "TRUE" = 2.2), guide = "none") +
  scale_alpha_identity() +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_void(base_size = 10) +
  theme(
    plot.background   = element_rect(fill = "#1a1a1a", color = NA),
    panel.background  = element_rect(fill = "#1a1a1a", color = NA),
    axis.text.x       = element_text(size = 8, color = "#aaaaaa"),
    axis.text.y       = element_text(size = 8, color = "#aaaaaa"),
    axis.title        = element_text(size = 9, color = "#aaaaaa"),
    axis.ticks        = element_line(color = "#444444"),
    panel.grid.major  = element_line(color = "#2e2e2e", linewidth = 0.3),
    legend.position   = "none",   # plotly colorbar replaces ggplot legend
    plot.margin       = margin(8, 8, 8, 8)
  ) +
  labs(x = "Poverty Rate", y = "Unemployment Rate", subtitle = NULL)

ggplotly(p_scatter, tooltip = "text") %>%
  layout(
    title  = list(
      text = "Poverty vs. Unemployment: All AZ Tracts",
      font = list(size = 12, color = "#cccccc")
    ),
    paper_bgcolor = "#1a1a1a",
    plot_bgcolor  = "#1a1a1a",
    xaxis = list(tickfont = list(color = "#aaaaaa"), titlefont = list(color = "#aaaaaa"),
                 gridcolor = "#2e2e2e", zerolinecolor = "#444444"),
    yaxis = list(tickfont = list(color = "#aaaaaa"), titlefont = list(color = "#aaaaaa"),
                 gridcolor = "#2e2e2e", zerolinecolor = "#444444"),
    showlegend = FALSE,
    coloraxis  = list(
      colorscale = list(
        list(0,   "#1696d2"),
        list(0.5, "#2d2d2d"),
        list(1,   "#db2b27")
      ),
      cmin = -3, cmax = 3,
      colorbar = list(showscale = FALSE)
    ),
    annotations = list(list(
      xref      = "paper", yref = "paper",
      x = 0.01, y = 0.99,
      xanchor = "left", yanchor = "top",
      align   = "left",
      showarrow = FALSE,
      bgcolor   = "rgba(15,15,15,0.82)",
      borderpad = 14,
      font      = list(color = "#aaaaaa", size = 13, family = "Lato"),
      text = paste0(
        "<b style='color:#ffffff;font-size:13px;letter-spacing:0.03em;'>Economic Hardship Index</b><br><br>",
        "<span style='color:#73bfe2;font-size:11px;font-weight:600;'>Better-off&nbsp;</span>",
        "<span style='color:#1696d2;font-size:19px;'>█</span>",
        "<span style='color:#4baecb;font-size:19px;'>█</span>",
        "<span style='color:#73bfe2;font-size:19px;'>█</span>",
        "<span style='color:#a8cfe0;font-size:19px;'>█</span>",
        "<span style='color:#d8d8d8;font-size:19px;'>█</span>",
        "<span style='color:#e8a090;font-size:19px;'>█</span>",
        "<span style='color:#ec6a5a;font-size:19px;'>█</span>",
        "<span style='color:#e03d30;font-size:19px;'>█</span>",
        "<span style='color:#db2b27;font-size:19px;'>█</span>",
        "<span style='color:#ec6a5a;font-size:11px;font-weight:600;'>&nbsp;Worse-off</span>"
      )
    )),
    margin = list(t = 55, b = 40, l = 50, r = 10)
  ) %>%
  config(displayModeBar = FALSE)
```

### LISA Map {width=33%}

```{r}
#| label: map-az-lisa
#| title: "Economic Hardship Clusters (LISA)"
# Re-simplify at 100m (vs the 500m used globally) to close crack artifacts at closer zoom
az_lisa_map_fine <- az_lisa_map %>%
  st_simplify(preserveTopology = TRUE, dTolerance = 100) %>%
  st_transform(4326)

maplibre(
  style  = carto_style("dark-matter"),
  bounds = c(-115.8, 31.3, -110.1, 37.0)
) |>
  add_fill_layer(
    id           = "az-lisa",
    source       = az_lisa_map_fine,
    fill_color   = match_expr(
      column  = "lisa_az",
      values  = c("HH",      "LL",      "HL",      "LH",      "NS"),
      stops   = c("#db2b27", "#1696d2", "#fdbf11", "#73bfe2", "#d2d2d2"),
      default = "#d2d2d2"
    ),
    fill_opacity = 0.85,
    tooltip      = "tooltip_text"
  ) |>
  add_line_layer(
    id = "maricopa-border-lisa", source = maricopa_boundary,
    line_color = "#fdbf11", line_width = 2.5
  ) |>
  add_line_layer(
    id           = "az-lisa-borders",
    source       = az_lisa_map_fine,
    line_color   = "#2a2a2a",
    line_width   = 0.4,
    line_opacity = 0.6
  ) |>
  htmlwidgets::onRender("
    function(el, x) {
      var entries = [
        { color: '#db2b27', label: 'Hot Spot (HH)' },
        { color: '#1696d2', label: 'Cold Spot (LL)' },
        { color: '#fdbf11', label: 'High-Low Outlier (HL)' },
        { color: '#73bfe2', label: 'Low-High Outlier (LH)' },
        { color: '#d2d2d2', label: 'Not Significant (NS)' }
      ];
      var leg = document.createElement('div');
      leg.style.cssText = 'position:absolute;bottom:28px;left:10px;background:rgba(15,15,15,0.82);color:#e8e8e8;font-size:9px;font-family:Lato,sans-serif;padding:6px 8px;border-radius:4px;line-height:1.5;min-width:170px;';
      var title = document.createElement('div');
      title.style.cssText = 'font-weight:700;font-size:9.5px;margin-bottom:5px;color:#ffffff;border-bottom:1px solid #444;padding-bottom:4px;';
      title.textContent = 'LISA Clusters (p < 0.05)';
      leg.appendChild(title);
      entries.forEach(function(e) {
        var row = document.createElement('div');
        row.style.cssText = 'display:flex;align-items:center;margin-bottom:2px;';
        var swatch = document.createElement('span');
        swatch.style.cssText = 'display:inline-block;width:9px;height:9px;border-radius:2px;margin-right:6px;flex-shrink:0;background:' + e.color + ';';
        var lbl = document.createElement('span');
        lbl.textContent = e.label;
        row.appendChild(swatch); row.appendChild(lbl);
        leg.appendChild(row);
      });
      var badge = document.createElement('div');
      badge.style.cssText = 'position:absolute;top:10px;right:12px;background:rgba(10,10,10,0.75);color:#fdbf11;font-size:11px;font-family:Lato,sans-serif;padding:5px 9px;border-radius:3px;border:1px solid #fdbf11;pointer-events:none;';
      badge.innerHTML = '\u25a0&nbsp;Maricopa County boundary';
      el.appendChild(leg);
      el.appendChild(badge);
    }
  ")
```

## Bottom Row {height=55%}

### County Table {width=67%}

```{r}
#| label: table-county-hetero
#| title: "Within-County Variation in Economic Hardship (EH) — All 15 Arizona Counties"
county_name_lookup <- county_2023 %>%
  st_drop_geometry() %>%
  filter(str_detect(NAME, "Arizona")) %>%
  mutate(
    county_geoid = GEOID,
    county_name  = str_remove(str_extract(NAME, "^[^,]+"), " County")
  ) %>%
  select(county_geoid, county_name)

# Build extra component summary columns dynamically
extra_summary_cols <- if (length(STUDENT_COMPONENTS) > 0) {
  setNames(
    lapply(STUDENT_COMPONENTS, function(comp)
      rlang::quo(scales::percent(mean(!!rlang::sym(comp), na.rm = TRUE), accuracy = 0.1))
    ),
    sapply(STUDENT_COMPONENTS, function(comp) EXTRA_VAR_MENU[[comp]]$label)
  )
} else list()

az_county_hetero_tbl <- az_tract_2023 %>%
  st_drop_geometry() %>%
  filter(!is.na(hardship_index)) %>%
  mutate(county_geoid = substr(GEOID, 1, 5)) %>%
  left_join(county_name_lookup, by = "county_geoid") %>%
  filter(!is.na(county_name)) %>%
  group_by(county_name) %>%
  summarise(
    `N Tracts`    = n(),
    `Pov. Rate`   = scales::percent(mean(poverty_rate,  na.rm = TRUE), accuracy = 0.1),
    `Unemp. Rate` = scales::percent(mean(unemp_rate,    na.rm = TRUE), accuracy = 0.1),
    `Med. Income` = scales::dollar(mean(median_income,  na.rm = TRUE), accuracy = 1),
    !!!extra_summary_cols,
    `EHI Mean`    = round(mean(hardship_index,          na.rm = TRUE), 3),
    .groups       = "drop"
  ) %>%
  arrange(`EHI Mean`) %>%
  rename(County = county_name)

# Dynamic header: colspan expands with extra components
n_comp_cols <- 3 + length(STUDENT_COMPONENTS)
extra_th    <- lapply(STUDENT_COMPONENTS, function(comp) htmltools::tags$th(EXTRA_VAR_MENU[[comp]]$label))

az_county_hetero_tbl %>%
  DT::datatable(
    options = list(
      pageLength = 8, dom = "t",
      scrollY = "280px", scrollX = TRUE,
      columnDefs = list(list(className = "dt-center", targets = 1:(n_comp_cols + 1)))
    ),
    rownames = FALSE,
    container = htmltools::withTags(
      table(
        class = "display",
        thead(
          tr(
            th(rowspan = 2, "County"),
            th(rowspan = 2, "N Tracts"),
            th(colspan = n_comp_cols,
               style = "text-align:center;border-bottom:2px solid #ddd;",
               "Index Components (County Avg.)"),
            th(rowspan = 2,
               style = "text-align:center;border-left:2px solid #ccc;padding-left:8px;",
               "EH Index")
          ),
          tr(
            th("Pov. Rate"), th("Unemp. Rate"), th("Med. Income"),
            extra_th
          )
        )
      )
    ),
    caption = htmltools::tags$caption(
      style = "font-size:0.71rem;color:#767676;text-align:left;caption-side:bottom;padding-top:4px;",
      paste0("Index composition: ", INDEX_LABEL, ". Sorted best-off to worst-off.")
    )
  ) %>%
  DT::formatStyle("County", target = "row",
    backgroundColor = DT::styleEqual(TARGET_COUNTY, "#fffde6"))
```

### AZ Alert Stats {width=33%}

#### Row

```{r}
#| label: vb-az-hh
#| content: valuebox
list(
  title = "High Economic Hardship Clusters (HH) — Hot Spot Tracts Statewide:",
  value = az_n_hh,
  icon  = "fire",
  color = "danger"
)
```

#### Row

```{r}
#| label: vb-az-ll
#| content: valuebox
list(
  title = "Low Economic Hardship Clusters (LL) — Cold Spot Tracts Statewide:",
  value = az_n_ll,
  icon  = "snow",
  color = "primary"
)
```



# Maricopa: Clusters & Trajectories

## Column {width=50%}


```{r}
#| label: sankey-ehi-quintiles
#| title: "Economic Hardship Mobility: Maricopa Tracts (2013 \u2192 2023)"
# Quintiles (5 bins) — cleaner than deciles for Sankey readability
ehi_flow <- change_df %>%
  st_drop_geometry() %>%
  select(GEOID, hardship_2013) %>%
  left_join(
    local_2023 %>% st_drop_geometry() %>% select(GEOID, hardship_2023 = hardship_index),
    by = "GEOID"
  ) %>%
  filter(!is.na(hardship_2013), !is.na(hardship_2023)) %>%
  mutate(
    q_2013 = ntile(hardship_2013, 5),
    q_2023 = ntile(hardship_2023, 5)
  ) %>%
  group_by(q_2013, q_2023) %>%
  summarise(n_tracts = n(), .groups = "drop")

q_labels <- c(
  "Q1  Least Hardship",
  "Q2",
  "Q3  Middle",
  "Q4",
  "Q5  Most Hardship"
)

# Blue (low) -> grey (mid) -> red (high)
q_colors <- colorRampPalette(c("#1696d2", "#d2d2d2", "#db2b27"))(5)

node_labels <- c(paste0("2013  ", q_labels), paste0("2023  ", q_labels))
node_colors <- c(q_colors, q_colors)

src <- ehi_flow$q_2013 - 1L
tgt <- ehi_flow$q_2023 - 1L + 5L

link_colors <- sapply(ehi_flow$q_2013, function(q) {
  h <- q_colors[q]
  r <- strtoi(substr(h, 2, 3), 16L)
  g <- strtoi(substr(h, 4, 5), 16L)
  b <- strtoi(substr(h, 6, 7), 16L)
  sprintf("rgba(%d,%d,%d,0.38)", r, g, b)
})

plot_ly(
  type        = "sankey",
  orientation = "h",
  arrangement = "freeform",
  node = list(
    label     = node_labels,
    color     = node_colors,
    pad       = 24,
    thickness = 28,
    line      = list(color = "#cccccc", width = 0.4)
  ),
  link = list(
    source = src,
    target = tgt,
    value  = ehi_flow$n_tracts,
    color  = link_colors,
    label  = paste0(ehi_flow$n_tracts, " tracts | ",
                    "2013 Q", ehi_flow$q_2013,
                    " \u2192 2023 Q", ehi_flow$q_2023)
  )
) %>%
  layout(
    title = list(
      text = paste0(
        "Economic Hardship Quintile Mobility: Maricopa County Tracts (2013 \u2192 2023)",
        "<br><sup>Q1 = least hardship &nbsp;\u2502&nbsp; Q5 = most hardship",
        " &nbsp;\u2502&nbsp; Width = tract count &nbsp;\u2502&nbsp; Hover for details</sup>"
      ),
      font = list(size = 12, color = "#333333")
    ),
    font   = list(size = 11, color = "#333333"),
    margin = list(t = 70, b = 20, l = 20, r = 20)
  ) %>%
  config(displayModeBar = FALSE)
```


## Column {width=50%}


### Row {height=60%}

```{r}
#| label: stats-2023
#| include: false
# Compute 2013->2023 improvement stats dynamically from actual data
maricopa_change_23 <- change_df %>%
  st_drop_geometry() %>%
  select(GEOID, hardship_2013, trajectory) %>%
  left_join(
    local_2023 %>% st_drop_geometry() %>% select(GEOID, hardship_2023 = hardship_index),
    by = "GEOID"
  ) %>%
  filter(!is.na(hardship_2013), !is.na(hardship_2023)) %>%
  mutate(change_13_23 = hardship_2023 - hardship_2013)

pct_improved_23   <- round(100 * mean(maricopa_change_23$change_13_23 < -0.05, na.rm = TRUE), 1)
pct_worsened_23   <- round(100 * mean(maricopa_change_23$change_13_23 >  0.05, na.rm = TRUE), 1)
pct_persist_hh_23 <- round(100 * mean(maricopa_change_23$trajectory == "Persistent HH", na.rm = TRUE), 1)
pct_emerg_hh_23   <- round(100 * mean(maricopa_change_23$trajectory == "Emerging HH",   na.rm = TRUE), 1)
```

```{r}
#| label: map-trajectory-gl
#| title: "Neighborhood Hardship Trajectories (2013\u20132019)"
traj_map <- change_df %>%
  mutate(
    trajectory   = as.character(trajectory),
    traj_tooltip = paste0(
      str_extract(NAME, "Census Tract [\\d\\.]+"), "\n",
      trajectory, "\n",
      "2013 EHI: ", round(hardship_2013, 3), "\n",
      "2019 EHI: ", round(hardship_2019, 3)
    )
  ) %>%
  st_simplify(preserveTopology = TRUE, dTolerance = 200) %>%
  st_transform(4326)

maplibre(
  style  = carto_style("dark-matter"),
  bounds = c(-113.3, 32.9, -111.0, 34.05)
) |>
  add_fill_layer(
    id           = "traj-fill",
    source       = traj_map,
    fill_color   = match_expr(
      column  = "trajectory",
      values  = c("Persistent HH", "Emerging HH",  "Dissolving HH",
                  "Persistent LL", "Emerging LL",   "HL Outlier",   "Stable NS"),
      stops   = c("#db2b27",       "#fdbf11",       "#fdd870",
                  "#1696d2",       "#73bfe2",        "#ec008b",      "#1a1a1a"),
      default = "#d2d2d2"
    ),
    fill_opacity = 0.85,
    tooltip      = "traj_tooltip"
  ) |>
  add_line_layer(
    id = "traj-border", source = traj_map,
    line_color = "#ffffff", line_width = 0.3, line_opacity = 0.3
  ) |>
  add_line_layer(
    id           = "maricopa-traj-boundary",
    source       = maricopa_boundary,
    line_color   = "#aaaaaa",
    line_width   = 1.8,
    line_opacity = 0.7
  ) |>
  htmlwidgets::onRender("
    function(el, x) {
      // Vertical custom legend — replaces mapgl default
      var entries = [
        { color: '#db2b27', label: 'Persistent Hot Spot (HH \u2192 HH)' },
        { color: '#fdbf11', label: 'Emerging Hot Spot (NS \u2192 HH)' },
        { color: '#fdd870', label: 'Dissolving Hot Spot (HH \u2192 NS)' },
        { color: '#1696d2', label: 'Persistent Cold Spot (LL \u2192 LL)' },
        { color: '#73bfe2', label: 'Emerging Cold Spot (NS \u2192 LL)' },
        { color: '#ec008b', label: 'Spatial Outlier (HL / LH)' }
      ];
      var leg = document.createElement('div');
      leg.style.cssText = 'position:absolute;top:12px;left:12px;background:rgba(15,15,15,0.82);color:#e8e8e8;font-size:11px;font-family:Lato,sans-serif;padding:10px 13px;border-radius:5px;line-height:1.7;min-width:210px;';
      var title = document.createElement('div');
      title.style.cssText = 'font-weight:700;font-size:11.5px;margin-bottom:7px;color:#ffffff;border-bottom:1px solid #444;padding-bottom:5px;';
      title.textContent = 'Cluster Trajectory (2013\u21922019)';
      leg.appendChild(title);
      entries.forEach(function(e) {
        var row = document.createElement('div');
        row.style.cssText = 'display:flex;align-items:center;margin-bottom:3px;';
        var swatch = document.createElement('span');
        swatch.style.cssText = 'display:inline-block;width:13px;height:13px;border-radius:2px;margin-right:8px;flex-shrink:0;background:' + e.color + ';';
        var lbl = document.createElement('span');
        lbl.textContent = e.label;
        row.appendChild(swatch); row.appendChild(lbl);
        leg.appendChild(row);
      });
      el.appendChild(leg);
    }
  ")
```

### Row {height=40%}

```{r}
#| label: vb-improved
#| content: valuebox
list(title = "Tracts Improved (2013\u21922023)",
     value = paste0(pct_improved_23, "%"),
     icon  = "arrow-up-circle", color = "primary")
```

```{r}
#| label: vb-worsened
#| content: valuebox
list(title = "Tracts Worsened (2013\u21922023)",
     value = paste0(pct_worsened_23, "%"),
     icon  = "arrow-down-circle", color = "warning")
```

```{r}
#| label: vb-persistent
#| content: valuebox
list(title = "Persistently High Hardship (2013\u20132023)",
     value = paste0(pct_persist_hh_23, "%"),
     icon  = "exclamation-circle", color = "danger")
```

```{r}
#| label: vb-emerging
#| content: valuebox
list(title = "Emerging Hot Spots (2013\u20132023)",
     value = paste0(pct_emerg_hh_23, "%"),
     icon  = "arrow-up-right-circle", color = "warning")
```

# Policy Implications

## Column {width=50%}

::: {.card style="border-left: 5px solid #db2b27; margin-bottom: 12px;"}
### 🔴 Areas of Persistent Concern

**`r n_persist_hh` Persistent Hot Spot Tracts (`r pct_persist_hh_23`% of Maricopa tracts)**

These census tracts had statistically significant high-hardship clustering in **both** 2013 and 2019. The pattern is not random: a Global Moran's I of `r round(moran_global$estimate[1], 3)` confirms that hardship is spatially concentrated, not scattered. Tracts in this category share infrastructure deficits, limited employment access, and concentrated poverty that reinforce one another across neighborhood boundaries.

**Implication:** Individual-level interventions alone are unlikely to move the needle. Place-based, multi-sector investment is required.

*[TODO: Identify the specific geographic corridor in your county where these tracts are concentrated.]*
:::

::: {.card style="border-left: 5px solid #fdbf11; margin-bottom: 12px;"}
### 🟡 Early Warning Signals

**`r n_emerg_hh` Emerging Hot Spot Tracts (`r pct_emerg_hh_23`% of Maricopa tracts)**

These tracts were **not** significant hardship clusters in 2013 but became statistically significant by 2019, representing the spatial expansion of hardship beyond historically distressed cores. This is an early warning signal that hardship is spreading, not contained.

**Displacement paradox:** Some "improving" tracts nearby may be gentrifying, pushing lower-income households outward into these emerging clusters. Declining hardship scores do not necessarily mean existing residents are better off.

*[TODO: Identify where emerging hot spots are forming in your county and what may be driving displacement.]*
:::

::: {.card style="border-left: 5px solid #55b748;"}
### 🟢 Signs of Progress Read With Caution

**`r pct_improved_23`% of Maricopa tracts showed EHI improvement (2013→2023)**

The majority of tracts improved over the decade-long window. However, aggregate improvement masks significant variation: `r pct_worsened_23`% of tracts worsened over the same period. The data cannot distinguish genuine economic uplift from population turnover: a tract with a declining hardship index may simply have replaced lower-income residents with higher-income newcomers.

**Data limitation:** Before drawing conclusions from improving scores, ground-truth verification through community engagement and displacement tracking is essential.

*[TODO: Customize with your county's specific context.]*
:::

## Column {width=50%}

::: {.card}
### Recommendation 1
**Target:** `r n_persist_hh` Persistent Hot Spot tracts (`r pct_persist_hh_23`% of all Maricopa tracts): spatially concentrated, entrenched hardship confirmed by a Global Moran's I of `r round(moran_global$estimate[1], 3)`.

**TODO:** Write your first recommendation. Name the specific geographic corridor, cite the hardship index values, and propose a concrete place-based intervention with a named responsible entity.

*[Your answer here minimum 3 sentences]*
:::

::: {.card}
### Recommendation 2
**Target:** `r n_emerg_hh` Emerging Hot Spot tracts (`r pct_emerg_hh_23`% of all Maricopa tracts), new high-hardship clusters not present in 2013, signaling spatial expansion.

**TODO:** Reference the displacement paradox, identify where these clusters are forming, and propose an early-intervention or monitoring strategy.

*[Your answer here minimum 3 sentences]*
:::

::: {.card}
### Recommendation 3
**Evidence base:** `r pct_improved_23`% of tracts improved (2013→2023) but `r pct_worsened_23`% worsened; Moran's I = `r round(moran_global$estimate[1], 3)` confirms strong spatial clustering persists.

**TODO:** Using the trajectory map and mobility Sankey, make a forward-looking data-monitoring or cross-sector coordination argument.

*[Your answer here minimum 3 sentences]*
:::

# Index Sensitivity

## Column {width=100%}

::: {.card style="border-left: 5px solid #9370db;"}
### 🔬 Index Sensitivity Reflection

The baseline EHI consists of **3 measures**: Poverty + Unemployment + Income (inv.)

**Current index:** `r INDEX_N`-component EHI: `r INDEX_LABEL`

**After adding your extra component(s), answer the following (minimum 2 sentences each):**

---

**Q1: What changed spatially?**  
Compare Hot Spot tract counts and cluster map patterns between your expanded index and the 3-component baseline. Did adding `r if(length(STUDENT_COMPONENTS)>0) paste(sapply(STUDENT_COMPONENTS, function(c) EXTRA_VAR_MENU[[c]]$label), collapse=" + ") else "[your chosen component]"` shift which tracts or corridors are flagged?

*[TODO: Describe what changed in the spatial cluster map and Hot Spot tract counts after adding your component. Be specific cite numbers and identify geographic areas.]*

**Q2: What stayed the same?**  
Which Persistent Hot Spot areas appear robustly across index specifications? What does consistency across different index compositions tell us about the reliability of hardship diagnoses in those tracts?

*[TODO: Identify which Persistent Hot Spot corridors remained flagged regardless of index composition. Explain what this robustness tells us about entrenched hardship in those areas.]*

**Q3: Policy implications of index choice**  
If a policymaker targeted place-based investments using the baseline index versus your expanded index, would resource allocation differ? Name specific tracts or geographic corridors and argue which composition better captures the full burden of economic hardship for policy purposes.

*[TODO: Compare resource targeting under each index. Argue which composition is more appropriate for policy, citing specific tracts or neighborhoods and the dimension of hardship your added component captures.]*
:::
