title: “lab7” author: “Joshua Hawk” date: “2026-03-26” output: html_document

As new people are moving to central Tennessee year by year means that there is a whole mixture of different languages and cultures coming to central Tennessee. Below is the code for the map and the graph that shows the difference in pupulation of people who speak languages other than english over a 5 year period (2019-2024).

# ----------------------------------------------------------
# Step 1: Install required packages (if missing)
# ----------------------------------------------------------

if (!require("tidyverse"))
  install.packages("tidyverse")
if (!require("tidycensus"))
  install.packages("tidycensus")
if (!require("sf"))
  install.packages("sf")
if (!require("leaflet"))
  install.packages("leaflet")
if (!require("htmlwidgets"))
  install.packages("htmlwidgets")
if (!require("plotly"))
  install.packages("plotly")

# ----------------------------------------------------------
# Step 2: Load libraries
# ----------------------------------------------------------

library(tidyverse)
library(tidycensus)
library(sf)
library(leaflet)
library(htmlwidgets)
library(plotly)

# ----------------------------------------------------------
# Step 3: Transmit Census API key
# ----------------------------------------------------------

census_api_key("494018a70f114d4f76b10537730ccc9c7dbfe36b", install = FALSE)

# ----------------------------------------------------------
# Step 4: Define analysis years, variables, labels, and codebooks
# ----------------------------------------------------------

year_old <- 2019  # Avoid overlapping ACS datasets
year_new <- 2024  # Avoid overlapping ACS datasets

# ---- Load codebooks for both years (for inspection / validation) ----

# Detailed tables
Old_DetailedTables <- load_variables(year_old, "acs5", cache = TRUE)
New_DetailedTables <- load_variables(year_new, "acs5", cache = TRUE)

# Subject tables
Old_SubjectTables  <- load_variables(year_old, "acs5/subject", cache = TRUE)
New_SubjectTables  <- load_variables(year_new, "acs5/subject", cache = TRUE)

# Profile tables
Old_ProfileTables  <- load_variables(year_old, "acs5/profile", cache = TRUE)
New_ProfileTables  <- load_variables(year_new, "acs5/profile", cache = TRUE)

# ---- Enter variable name(s); names often change year to year ----
Variable_old <- c(Estimate = "DP02_0117P")
Variable_new <- c(Estimate = "DP02_0117P")

# Human-readable label for variable
target_label <- "Total Population" # <<< Edit as needed

# ----------------------------------------------------------
# Step 5: Fetch ACS data – newer year
# ----------------------------------------------------------

mydata_new <- get_acs(
  geography = "county subdivision",
  state = "TN",
  variables = Variable_new,
  year = year_new,
  survey = "acs5",
  output = "wide",
  geometry = TRUE
) %>%
  rename(Area = NAME) %>%
  transmute(
    GEOID,
    Area,
    geometry,
    Estimate_new = EstimateE,
    MoE_new      = EstimateM
  )

# ----------------------------------------------------------
# Step 6: Fetch ACS data – older year
# ----------------------------------------------------------

mydata_old <- get_acs(
  geography = "county subdivision",
  state = "TN",
  variables = Variable_old,
  year = year_old,
  survey = "acs5",
  output = "wide",
  geometry = FALSE
) %>%
  rename(Area = NAME) %>%
  transmute(
    GEOID,
    Area,
    Estimate_old = EstimateE,
    MoE_old      = EstimateM
  )

# ----------------------------------------------------------
# Step 7: Join datasets
# ----------------------------------------------------------

mydata_joined <- mydata_new %>%
  left_join(mydata_old, by = c("GEOID", "Area"))

# ----------------------------------------------------------
# Step 8: Flexible filtering by Area
# ----------------------------------------------------------

search_terms <- c("Rutherford", "Davidson", "Williamson")
# Use c("ALL") to include all areas

if (length(search_terms) == 1 && toupper(search_terms) == "ALL") {
  filtereddata <- mydata_joined
} else {
  or_pattern <- paste(
    stringr::str_replace_all(search_terms, "([\\W])", "\\\\\\1"),
    collapse = "|"
  )
  
  filtereddata <- mydata_joined %>%
    filter(
      stringr::str_detect(
        Area %>% replace_na(""),
        regex(or_pattern, ignore_case = TRUE)
      )
    )
}

# ----------------------------------------------------------
# Step 9: Compute change and statistical significance
# ----------------------------------------------------------

filtereddata <- filtereddata %>%
  mutate(
    Diff     = Estimate_new - Estimate_old,
    Diff_MOE = tidycensus::moe_sum(MoE_new, MoE_old),
    Sig_90   = tidycensus::significance(
      Estimate_new, MoE_new,
      Estimate_old, MoE_old
    ),
    SigDirection = case_when(
      Sig_90 & Diff > 0 ~ "Significant increase",
      Sig_90 & Diff < 0 ~ "Significant decrease",
      TRUE             ~ "No statistically significant change"
    )
  )

# ----------------------------------------------------------
# Step 10: Prepare spatial data
# ----------------------------------------------------------

mapdata <- filtereddata %>%
  st_as_sf() %>%
  st_transform(4326)

# ----------------------------------------------------------
# Step 11: Ranked change dot plot (significance encoded by color)
# ----------------------------------------------------------

plotdf <- filtereddata %>%
  st_drop_geometry() %>%
  mutate(
    y_ordered = reorder(Area, Diff),
    sig_class = case_when(
      SigDirection == "Significant increase" ~ "Increase (significant)",
      SigDirection == "Significant decrease" ~ "Decrease (significant)",
      TRUE ~ "Not statistically significant"
    ),
    hover_text = paste0(
      "Area: ", Area,
      "<br>Change (", year_old, " to ", year_new, "): ", scales::comma(Diff),
      "<br>Significance: ", sig_class,
      "<br>", year_old, ": ", scales::comma(Estimate_old),
      "<br>", year_new, ": ", scales::comma(Estimate_new)
    )
  )

mygraph <- plotly::plot_ly(
  data = plotdf,
  x = ~Diff,
  y = ~as.character(y_ordered),
  type = "scatter",
  mode = "markers",
  color = ~sig_class,
  colors = c(
    "Increase (significant)" = "#b2182b",
    "Decrease (significant)" = "#2166ac",
    "Not statistically significant" = "#bdbdbd"
  ),
  marker = list(size = 9),
  text = ~hover_text,
  hovertemplate = "%{text}<extra></extra>"
) %>%
  layout(
    title = list(
      text = paste0(
        "Estimated Change in ", target_label,
        " (", year_old, " to ", year_new, ")",
        "<br><sup>County subdivisions. Color indicates statistical significance (90%).</sup>"
      )
    ),
    xaxis = list(title = "Estimated change", tickformat = ",.0f"),
    yaxis = list(title = "", automargin = TRUE)
  )

mygraph

# ----------------------------------------------------------
# Step 12: Leaflet popup content (no MOE for change)
# ----------------------------------------------------------

mapdata$popup <- paste0(
  "<strong>", mapdata$Area, "</strong><br/>",
  "<hr>",
  year_old, ": ", format(mapdata$Estimate_old, big.mark = ","),
  " (plus/minus ", format(mapdata$MoE_old, big.mark = ","), ")<br/>",
  year_new, ": ", format(mapdata$Estimate_new, big.mark = ","),
  " (plus/minus ", format(mapdata$MoE_new, big.mark = ","), ")<br/>",
  "<b>Change:</b> ", format(mapdata$Diff, big.mark = ","), "<br/>",
  "<em>", mapdata$SigDirection, " (90%)</em>"
)

# ----------------------------------------------------------
# Step 13: Leaflet map of change
# ----------------------------------------------------------

max_abs <- max(abs(mapdata$Diff), na.rm = TRUE)
bins <- pretty(c(-max_abs, max_abs), n = 6)

pal <- colorBin(
  palette = "RdBu",
  domain  = mapdata$Diff,
  bins    = bins,
  reverse = TRUE
)

DivisionMap <- leaflet(mapdata) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor   = ~pal(Diff),
    fillOpacity = 0.6,
    color       = "black",
    weight      = 1,
    popup       = ~popup
  ) %>%
  addLegend(
    pal    = pal,
    values = ~Diff,
    title  = paste0(
      "Change (", year_old, " to ", year_new, ")",
      "<br/><small>", target_label, "</small>"
    ),
    labFormat = labelFormat(big.mark = ",")
  )

DivisionMap

# ----------------------------------------------------------
# Step 14: Export outputs
# ----------------------------------------------------------

saveWidget(as_widget(mygraph), "ACSGraph_Diff.html", selfcontained = TRUE)
saveWidget(DivisionMap, "ACSMap_Diff.html", selfcontained = TRUE)