library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'ggplot2' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'stringr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(usmap)
library(viridis)
## Warning: package 'viridis' was built under R version 4.3.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.3.2
# Load data
df <- read.csv("C:/Users/mferdo2/OneDrive - Louisiana State University/Finance_PhD/Real_Estate_project/02_data/processed/broadband_annual_2018_2024.csv")

# Calculate broadband penetration
df <- df %>%
  mutate(
    fips = sprintf("%05d", county_fips),
    bbd_penetration = (consumer_bbd / Housing_Units) * 100
  )
# ===== TOP/BOTTOM 5% =====

# Filter 2024 data and calculate percentiles
df_2024 <- df %>%
  filter(year == 2024, !is.na(bbd_penetration))

# Get thresholds
top_5_threshold <- quantile(df_2024$bbd_penetration, 0.95, na.rm = TRUE)
bottom_5_threshold <- quantile(df_2024$bbd_penetration, 0.05, na.rm = TRUE)

# Create category variable
df_2024 <- df_2024 %>%
  mutate(
    category = case_when(
      bbd_penetration >= top_5_threshold ~ "Top 5%",
      bbd_penetration <= bottom_5_threshold ~ "Bottom 5%",
      TRUE ~ "Middle 90%"
    )
  )

# Summary stats
cat("Top 5% threshold:", round(top_5_threshold, 1), "\n")
## Top 5% threshold: 0.1
cat("Bottom 5% threshold:", round(bottom_5_threshold, 1), "\n")
## Bottom 5% threshold: 0
cat("Counties in top 5%:", sum(df_2024$category == "Top 5%"), "\n")
## Counties in top 5%: 159
cat("Counties in bottom 5%:", sum(df_2024$category == "Bottom 5%"), "\n")
## Counties in bottom 5%: 159
# Map - Top/Bottom 5%
plot_usmap(
  data = df_2024, 
  values = "category",
  regions = "counties",
  color = "white",
  size = 0.1
) +
  scale_fill_manual(
    name = "Broadband Penetration",
    values = c(
      "Top 5%" = "#2a9d8f",      # teal
      "Bottom 5%" = "#e63946",   # red
      "Middle 90%" = "gray85"
    ),
    breaks = c("Top 5%", "Middle 90%", "Bottom 5%")
  ) +
  labs(
    title = "Top 5% and Bottom 5% Counties by Broadband Penetration (2024)",
    subtitle = "Consumer broadband connections per 100 housing units",
    caption = paste0("Top 5% threshold: ", round(top_5_threshold, 1), 
                    " | Bottom 5% threshold: ", round(bottom_5_threshold, 1))
  ) +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    legend.position = "right",
    legend.text = element_text(size = 10)
  )

ggsave("broadband_top_bottom_5pct.png", width = 12, height = 8, dpi = 300)
# =====TOP/BOTTOM 1% =====

# Get thresholds
top_1_threshold <- quantile(df_2024$bbd_penetration, 0.99, na.rm = TRUE)
bottom_1_threshold <- quantile(df_2024$bbd_penetration, 0.01, na.rm = TRUE)

# Create category variable
df_2024 <- df_2024 %>%
  mutate(
    category_1pct = case_when(
      bbd_penetration >= top_1_threshold ~ "Top 1%",
      bbd_penetration <= bottom_1_threshold ~ "Bottom 1%",
      TRUE ~ "Middle 98%"
    )
  )

# Summary stats
cat("\nTop 1% threshold:", round(top_1_threshold, 1), "\n")
## 
## Top 1% threshold: 0.1
cat("Bottom 1% threshold:", round(bottom_1_threshold, 1), "\n")
## Bottom 1% threshold: 0
cat("Counties in top 1%:", sum(df_2024$category_1pct == "Top 1%"), "\n")
## Counties in top 1%: 32
cat("Counties in bottom 1%:", sum(df_2024$category_1pct == "Bottom 1%"), "\n")
## Counties in bottom 1%: 49
# Map - Top/Bottom 1%
plot_usmap(
  data = df_2024, 
  values = "category_1pct",
  regions = "counties",
  color = "white",
  size = 0.1
) +
  scale_fill_manual(
    name = "Broadband Penetration",
    values = c(
      "Top 1%" = "#2a9d8f",      # teal
      "Bottom 1%" = "#e63946",   # red
      "Middle 98%" = "gray85"
    ),
    breaks = c("Top 1%", "Middle 98%", "Bottom 1%")
  ) +
  labs(
    title = "Top 1% and Bottom 1% Counties by Broadband Penetration (2024)",
    subtitle = "Consumer broadband connections per 100 housing units",
    caption = paste0("Top 1% threshold: ", round(top_1_threshold, 1), 
                    " | Bottom 1% threshold: ", round(bottom_1_threshold, 1))
  ) +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    legend.position = "right",
    legend.text = element_text(size = 10)
  )

ggsave("broadband_top_bottom_1pct.png", width = 12, height = 8, dpi = 300)
# ===== EXPORT COUNTY LISTS =====

# Top 5% counties
top_5pct <- df_2024 %>%
  filter(category == "Top 5%") %>%
  arrange(desc(bbd_penetration)) %>%
  select(fips, countyname, statename, bbd_penetration, Housing_Units, consumer_bbd)

write.csv(top_5pct, "top_5pct_counties.csv", row.names = FALSE)

# Bottom 5% counties
bottom_5pct <- df_2024 %>%
  filter(category == "Bottom 5%") %>%
  arrange(bbd_penetration) %>%
  select(fips, countyname, statename, bbd_penetration, Housing_Units, consumer_bbd)

write.csv(bottom_5pct, "bottom_5pct_counties.csv", row.names = FALSE)

# Top 1% counties
top_1pct <- df_2024 %>%
  filter(category_1pct == "Top 1%") %>%
  arrange(desc(bbd_penetration)) %>%
  select(fips, countyname, statename, bbd_penetration, Housing_Units, consumer_bbd)

write.csv(top_1pct, "top_1pct_counties.csv", row.names = FALSE)

# Bottom 1% counties
bottom_1pct <- df_2024 %>%
  filter(category_1pct == "Bottom 1%") %>%
  arrange(bbd_penetration) %>%
  select(fips, countyname, statename, bbd_penetration, Housing_Units, consumer_bbd)

write.csv(bottom_1pct, "bottom_1pct_counties.csv", row.names = FALSE)

cat("\nFiles saved:\n")
## 
## Files saved:
cat("- broadband_top_bottom_5pct.png\n")
## - broadband_top_bottom_5pct.png
cat("- broadband_top_bottom_1pct.png\n")
## - broadband_top_bottom_1pct.png
cat("- top_5pct_counties.csv\n")
## - top_5pct_counties.csv
cat("- bottom_5pct_counties.csv\n")
## - bottom_5pct_counties.csv
cat("- top_1pct_counties.csv\n")
## - top_1pct_counties.csv
cat("- bottom_1pct_counties.csv\n")
## - bottom_1pct_counties.csv