In ctesta01/ggtriangles: Implements geom_triangles

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 18, 
  fig.height = 14
)

library(ggtriangles)
library(urbnmapr)


# download county and state maps
counties <- get_urbn_map('counties', sf=T)
states <- get_urbn_map('states', sf=T)

darkgreen <- colorspace::lighten('#4e7c71')
lightgreen <- colorspace::lighten("#bacac9", amount = 0.5)

Now do it with NYT data

nyt2020 <- readr::read_csv("https://github.com/nytimes/covid-19-data/raw/master/rolling-averages/us-counties-2020.csv")
nyt2021 <- readr::read_csv("https://github.com/nytimes/covid-19-data/raw/master/rolling-averages/us-counties-2021.csv") 
nyt_anomalies <- readr::read_csv("https://raw.githubusercontent.com/nytimes/covid-19-data/master/rolling-averages/anomalies.csv")

nyt2020avg <- nyt2020 %>% 
  bind_rows(nyt2021 %>% 
  filter(date <= lubridate::ymd("2021-04-19"))) %>% 
  left_join(
    nyt_anomalies %>% 
      filter(type %in% c('deaths', 'both')),
    by = c('geoid' = 'geoid', 'date' = 'date', 'county' = 'county', 'state' = 'state')
  ) %>% 
  filter(is.na(omit_from_rolling_average_on_subgeographies) | omit_from_rolling_average_on_subgeographies != 'yes') %>% 
  filter(date > lubridate::ymd("2020-04-19")) %>% 
  group_by(geoid, county, state) %>% 
  summarize(avg_death_rate = mean(deaths_avg_per_100k, na.rm=T))

nyt2021avg <- nyt2021 %>% group_by(geoid, county, state) %>% 
  filter(date > lubridate::ymd("2021-04-19")) %>% 
  filter(date <= lubridate::ymd("2021-12-28")) %>% 
 left_join(
    nyt_anomalies %>% 
      filter(type %in% c('deaths', 'both')),
    by = c('geoid' = 'geoid', 'date' = 'date', 'county' = 'county', 'state' = 'state')
  ) %>% 
 filter(is.na(omit_from_rolling_average_on_subgeographies) | omit_from_rolling_average_on_subgeographies != 'yes') %>% 
  summarize(avg_death_rate = mean(deaths_avg_per_100k, na.rm=T))


library(magrittr)
nyt2020avg %<>% rename(avg2020 = avg_death_rate)
nyt2021avg %<>% rename(avg2021 = avg_death_rate)

avg_death_rates <- nyt2020avg %<>% left_join(nyt2021avg)
avg_death_rates %<>% mutate(change = avg2021 - avg2020)
avg_death_rates %<>% mutate(fips = stringr::str_extract(geoid, "[0-9]+"))

counties <- urbnmapr::get_urbn_map('counties', sf = TRUE)
counties %<>% left_join(avg_death_rates %>% select(fips, change), by = c('county_fips' = 'fips'))

state_abbreviations <- 
tibble::tribble(
         ~state_name, ~state_abbv, ~state_abbv_med,
         "Alabama",                 "AL",                 "Ala.",
          "Alaska",                 "AK",               "Alaska",
         "Arizona",                 "AZ",                "Ariz.",
        "Arkansas",                 "AR",                 "Ark.",
      "California",                 "CA",               "Calif.",
        "Colorado",                 "CO",               "Color.",
     "Connecticut",                 "CT",                "Conn.",
        "Delaware",                 "DE",                 "Del.",
         "Florida",                 "FL",                 "Fla.",
         "Georgia",                 "GA",                  "Ga.",
          "Hawaii",                 "HI",               "Hawaii",
           "Idaho",                 "ID",                "Idaho",
        "Illinois",                 "IL",                 "Ill.",
         "Indiana",                 "IN",                 "Ind.",
            "Iowa",                 "IA",                 "Iowa",
          "Kansas",                 "KS",                 "Kan.",
        "Kentucky",                 "KY",                  "Ky.",
       "Louisiana",                 "LA",                  "La.",
           "Maine",                 "ME",                "Maine",
        "Maryland",                 "MD",                  "Md.",
   "Massachusetts",                 "MA",                "Mass.",
        "Michigan",                 "MI",                "Mich.",
       "Minnesota",                 "MN",                "Minn.",
     "Mississippi",                 "MS",                "Miss.",
        "Missouri",                 "MO",                  "Mo.",
         "Montana",                 "MT",                "Mont.",
        "Nebraska",                 "NE",                 "Neb.",
          "Nevada",                 "NV",                 "Nev.",
   "New Hampshire",                 "NH",                 "N.H.",
      "New Jersey",                 "NJ",                 "N.J.",
      "New Mexico",                 "NM",                 "N.M.",
        "New York",                 "NY",                 "N.Y.",
  "North Carolina",                 "NC",                 "N.C.",
    "North Dakota",                 "ND",                 "N.D.",
            "Ohio",                 "OH",                 "Ohio",
        "Oklahoma",                 "OK",                "Okla.",
          "Oregon",                 "OR",                 "Ore.",
    "Pennsylvania",                 "PA",                  "Pa.",
    "Rhode Island",                 "RI",                 "R.I.",
  "South Carolina",                 "SC",                 "S.C.",
    "South Dakota",                 "SD",               "S.Dak.",
       "Tennessee",                 "TN",                "Tenn.",
           "Texas",                 "TX",                 "Tex.",
            "Utah",                 "UT",                 "Utah",
         "Vermont",                 "VT",                 "V.T.",
        "Virginia",                 "VA",                  "Va.",
      "Washington",                 "WA",                "Wash.",
   "West Virginia",                 "WV",                "W.Va.",
       "Wisconsin",                 "WI",                 "Wis.",
         "Wyoming",                 "WY",                 "Wyo."
  )


state_centers <- sf::st_centroid(states)
state_centers_xy <- sf::st_coordinates(state_centers)
state_centers %<>% bind_cols(state_centers_xy %>% as.data.frame())
state_centers %<>% left_join(state_abbreviations)

lower_threshold <- 0.01
upper_threshold <- -0.01

ggplot() +
  geom_sf(data = counties,
          mapping = aes(
          fill = factor(change > lower_threshold),
          size = factor(change > lower_threshold) %>% forcats::fct_explicit_na(na_level = "FALSE"),
          color = factor(change > lower_threshold)
          )) +
  geom_sf(data = states, aes(), size = 0.5, color = 'black', fill = NA) +
  geom_triangles(
    data =
      bind_cols(
        counties %>% filter(change > lower_threshold),
        counties %>% filter(change > lower_threshold) %>%
          sf::st_centroid() %>%
          sf::st_coordinates(.) %>%
          as.data.frame()),
    mapping = aes(x = X, y = Y, triangle_height = change),
    triangle_width = .075, lineend = 'butt',
    fill = '#c0392b', color = '#c0392b', size = 0) +
  geom_text(
    data = state_centers, 
    mapping = aes(x = X, y = Y, label = state_abbv_med),
    hjust = 0.5, vjust = 0, check_overlap = F, size = 2.5
  ) + 
  scale_fill_manual(values = c(`TRUE` = 'lightpink', `FALSE` = 'white')) +
  scale_size_manual(values = c(`TRUE` = .05, `FALSE` = 0)) +
  scale_color_manual(values = c(`TRUE` = "#c0392b", `FALSE` = 'white')) +
  scale_triangle_height(
    range = c(0.1, .75),
    n.breaks = 2,
    limits = c(0, max(abs(range(counties$change, na.rm=T)))),
    breaks = c(.1, 2.25),
    labels = c('smaller\nincrease', 'larger\nincrease')
    ) + 
  ggtitle("Counties where death rates were higher after vaccines were widely available") + 
  theme_void() + 
  labs(triangle_height = '') + 
  guides(triangle_height = guide_legend(order=1, label.position = 'bottom'), color = 'none', fill = 'none', size = 'none') + 
  theme(legend.position = 'bottom', plot.title = element_text(hjust = 0.5, face = 'bold'))


# library(here)
# ggsave(here("img/nyt_increased_counties.png"), width = 10, height = 8)
# ggsave(here("img/nyt_increased_counties.jpeg"), width = 10, height = 8)

ggplot() +
  geom_sf(data = counties,
          mapping = aes(
          fill = factor(change < upper_threshold),
          size = factor(change < upper_threshold) %>% forcats::fct_explicit_na(na_level = "FALSE"),
          color = factor(change < upper_threshold)
          )) +
  geom_sf(data = states, aes(), size = 0.5, color = 'black', fill = NA) +
  geom_triangles(
    data =
      bind_cols(
        counties %>% filter(change < upper_threshold),
        counties %>% filter(change < upper_threshold) %>%
          sf::st_centroid() %>%
          sf::st_coordinates(.) %>%
          as.data.frame()),
    mapping = aes(x = X, y = Y, triangle_height = change),
    triangle_width = .075, lineend = 'butt',
    fill = darkgreen, color = darkgreen,
    legend_y_offset = 0.5, size = 0) +
  geom_text(
    data = state_centers, 
    mapping = aes(x = X, y = Y, label = state_abbv_med),
    hjust = 0.5, vjust = 0, size = 2.5
  ) + 
  scale_triangle_height(
    range = c(-0.1,-.75),
    limits = c(-max(abs(range(counties$change, na.rm=T))), 0),
    n.breaks = 2,
    breaks = c(-.05,-2.25),
    labels = c('larger\ndecrease', 'smaller\ndecrease')
    ) + 
  scale_fill_manual(values = c(`TRUE` = lightgreen, `FALSE` = 'white')) +
  scale_size_manual(values = c(`TRUE` = .05, `FALSE` = 0)) +
  scale_color_manual(values = c(`TRUE` = darkgreen, `FALSE` = 'white')) +
  ggtitle("Counties where death rates were lower after vaccines were widely available") + 
  theme_void() + 
  labs(triangle_height = '') + 
  guides(triangle_height = guide_legend(order=1, reverse=TRUE, label.position = 'bottom'), color = 'none', fill = 'none', size = 'none') + 
  theme(legend.position = 'bottom',
        plot.title = element_text(hjust = 0.5, face = 'bold'))

# ggsave(here("img/nyt_decreased_counties.png"), width = 10, height = 8)
# ggsave(here("img/nyt_decreased_counties.jpeg"), width = 10, height = 8)