MetricT

New cases per 100k in US by state, facet

Jul 28th, 2020 (edited)
3,548
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ### Graph a geofacet map of new cases per 100k by state - /u/MetricT
  2.  
  3. library(tidyverse)
  4. library(tidycensus)
  5. library(cowplot)
  6. library(TTR)
  7. library(geofacet)
  8. library(broom)
  9.  
  10. ### Pull cases data from the NY Times repo below
  11. ###
  12. ### https://github.com/nytimes/covid-19-data
  13. ###
  14. spreadsheet <-
  15.   "../nytimes_covid-19-data/us-states.csv" %>%
  16.   read_csv(col_names = TRUE, col_types = "Dccdd")
  17.  
  18. ### How big should our SMA window be?
  19. SMA_n <- 7
  20.  
  21. ### Pull state population data from the Census
  22. pop2018 <-
  23.   get_acs(geography   = "state",
  24.           variables   = c("B01003_001"),
  25.           year        = 2018,
  26.           geometry    = FALSE,
  27.           cache_table = TRUE) %>%
  28.   rename(POP2018 = estimate) %>%
  29.   select(NAME, POP2018) %>%
  30.   left_join(fips_codes %>%
  31.               select(state_name, state) %>%
  32.               unique() %>%
  33.               as_tibble(),
  34.             by = c("NAME" = "state_name")) %>%
  35.   rename(state_abbr = state,
  36.          state = NAME) %>%
  37.   select(state, state_abbr, POP2018)
  38.  
  39. ### Take the data and convert it to a 7-day SMA of new cases per capita
  40. data <-
  41.   spreadsheet %>%
  42.   select(date, state, cases) %>%
  43.   arrange(date) %>%
  44.   pivot_wider(id_cols = "date", names_from = "state", values_from = "cases") %>%
  45.   mutate(across(!starts_with("date"),
  46.                 .fns = list(new = ~ c(.[1], diff(.))),
  47.                 .names = "{fn}_{col}"
  48.   )) %>%
  49.   select(date, contains("new_"))  %>%
  50.   rename_at(vars(contains("new_")),
  51.             ~ str_replace(., "new_", "")) %>%
  52.   mutate(across(!starts_with("date"),
  53.                 .fns = list(sma = ~ SMA(., n = SMA_n)),
  54.                 .names = "{fn}_{col}"
  55.   )) %>%
  56.   select(date, contains("sma_"))  %>%
  57.   rename_at(vars(contains("sma_")),
  58.             ~ str_replace(., "sma_", "")) %>%
  59.   pivot_longer(-date, names_to = "state", values_to = "new_cases") %>%
  60.   left_join(pop2018, by = c("state" = "state")) %>%
  61.   mutate(new_cases_percapita = 100000 * new_cases / POP2018) %>%
  62.   select(date, state, state_abbr, new_cases_percapita)
  63.  
  64.  
  65. ### Add dot to the end of the line for reference
  66. dots <-
  67.   data %>%
  68.   select(date, state, new_cases_percapita) %>%
  69.   pivot_wider(id_cols = "date",
  70.               names_from = "state",
  71.               values_from = "new_cases_percapita") %>%
  72.   tail(n = 1) %>%
  73.   pivot_longer(-date, names_to = "state", values_to = "new_cases_percapita") %>%
  74.   rename(dots = new_cases_percapita)
  75.  
  76. data <- data %>% left_join(dots, by = c("date" = "date", "state" = "state"))
  77.  
  78. ### Compute the growth rate over the last 2 weeks
  79. growth <-
  80.   data %>%
  81.   select(date, state, new_cases_percapita) %>%
  82.   filter(date >= as.Date(Sys.Date()) - 14) %>%
  83.   filter(!state %in% c("Virgin Islands",
  84.                        "Guam",
  85.                        "Northern Mariana Islands")) %>%
  86.   group_by(state) %>%
  87.   do(tidy(lm(new_cases_percapita ~ date, .))) %>%
  88.   filter(term == "date") %>%
  89.   select(state, estimate) %>%
  90.   rename(growth = estimate) %>%
  91.   arrange(growth)
  92.  
  93. #data <- data %>% left_join(dots, by = c("date" = "date", "state" = "state"))
  94.  
  95. this_title <-
  96.   paste("New Cases per 100k (", SMA_n, "-day moving average)", sep = "")
  97.  
  98. g_facet <-
  99.   ggplot(data = data) +
  100.   theme_bw() +
  101.   geom_line(aes(x = as.Date(date), y = new_cases_percapita)) +
  102.   geom_point(aes(x = as.Date(date), y = dots), color = "blue") +
  103.   facet_geo(~ state) +
  104.   labs(title = this_title,
  105.        caption = "", x = "", y = "")
  106. print(g_facet)
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×