· 6 years ago · Apr 21, 2020, 08:48 PM
1################################################################################
2### Render a map of Tennessee using TN Dept. of Health data on COVID-19
3### infections showing total cases, new cases, cases by age, and a map of
4### confirmed infected.
5###
6### TN Dept. of Health COVID-19 website:
7### https://www.tn.gov/health/cedep/ncov.html
8###
9### Note: You will need a Census API key to access their data. You
10### can get a key by going to this website:
11###
12### https://api.census.gov/data/key_signup.html
13###
14### Once you have the key, just add it below.
15################################################################################
16
17### Set the Census API key
18census_key <- "PUT_YOUR_CENSUS_API_KEY_HERE"
19
20### Let's start by cleaning the environment, it currently causes the graphs
21### some problems if the environment is already populated
22rm(list = ls())
23
24### Load the necessary libraries
25packages <- c("tidyverse", "dplyr", "readxl", "ggplot2", "sf", "rgeos", "TTR",
26 "scales", "cowplot", "viridis", "gridExtra", "tidycensus",
27 "googlesheets4")
28
29new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
30if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
31invisible(lapply(packages, "library", quietly = TRUE,
32 character.only = TRUE, warn.conflicts = FALSE))
33
34### Set the Census API key
35census_api_key(census_key)
36
37### Enable map caching
38options(tigris_use_cache = TRUE)
39
40
41################################################################################
42### Read data from my spreadsheet
43################################################################################
44### Load COVID data from Google Sheets
45gs4_deauth()
46covid_data <- "https://docs.google.com/spreadsheets/d/1Opqs3DIYTlPbPhqrUa66RfKBBDVo0u5sntFM6uYZDzk/edit?usp=sharing"
47infected_df <- read_sheet(covid_data, sheet = "Infected")
48deaths_df <- read_sheet(covid_data, sheet = "Deaths")
49hospital_df <- read_sheet(covid_data, sheet = "Hospitalizations")
50testres_df <- read_sheet(covid_data, sheet = "Test_Results")
51testrates_df <- read_sheet(covid_data, sheet = "Testing_Rates") %>%
52 mutate(County = gsub(" County", "", County))
53recover_df <- read_sheet(covid_data, sheet = "Recovered")
54age_df <- read_sheet(covid_data, sheet = "Age_Data")
55
56
57################################################################################
58### Subset the the data here (if desired)
59################################################################################
60### For the map, we only want the latest total for each county.
61non_county_fields <- c("Date", "Total", "Other_State_Country",
62 "Unknown", "ERROR!!!")
63
64### All counties
65all_counties <-
66 colnames(infected_df) %>%
67 enframe() %>%
68 filter(!colnames(infected_df) %in% non_county_fields) %>%
69 select("value") %>%
70 deframe()
71
72### Counties in the Nashville Metropolitan Statistical Area
73nashville_msa <-
74 c("Cannon", "Cheatham", "Davidson", "Dickson", "Macon", "Maury", "Robertson",
75 "Rutherford", "Smith", "Sumner", "Trousdale", "Williamson", "Wilson")
76
77middle_tn <- c(
78 "Bedford", "Cannon", "Cheatham", "Clay", "Coffee", "Davidson",
79 "DeKalb", "Dickson", "Fentress", "Franklin", "Giles", "Grundy",
80 "Hickman", "Houston", "Humphreys", "Jackson", "Lawrence", "Lewis",
81 "Lincoln", "Macon", "Marshall", "Maury", "Montgomery", "Moore",
82 "Overton", "Perry", "Pickett", "Putnam", "Robertson", "Rutherford",
83 "Sequatchie", "Smith", "Stewart", "Sumner", "Trousdale", "Van Buren",
84 "Warren", "Wayne", "White", "Williamson", "Wilson")
85
86
87################################################################################
88### Select the map you want
89################################################################################
90#graph_counties <- nashville_msa
91#map_counties <- nashville_msa
92#location <- "Nashville MSA"
93
94#graph_counties <- c("Davidson")
95#map_counties <- c("Davidson")
96#location <- "Davidson County"
97
98#graph_counties <- middle_tn
99#map_counties <- middle_tn
100#location <- "Middle TN"
101
102graph_counties <- c("Total")
103map_counties <- all_counties
104location <- "Tennessee"
105
106### Use TidyCensus to pull population data as well as map geometry
107### from the Census bureau
108county_acs <-
109 get_acs(geography = "county",
110 variables = c(POP2018 = "B01003_001"),
111 state = c("47"),
112 county = map_counties,
113 year = 2018,
114 geometry = TRUE,
115 cache_table = TRUE) %>%
116 rename(POP2018 = estimate)
117
118### Split the Census data apart into a map...
119tn_map_df <-
120 county_acs %>%
121 select(GEOID, geometry)
122
123### And a population dataframe
124tn_pop_df <-
125 county_acs %>%
126 st_set_geometry(NULL) %>%
127 mutate(NAME = str_replace(NAME, " County, Tennessee", "")) %>%
128 select(GEOID, NAME, POP2018) %>%
129 rename(County = NAME)
130
131
132################################################################################
133### A few visual configuration settings
134################################################################################
135map_palette <- c("#c4c4c4", viridis(n = 1000, direction = -1))
136map_fontsize <- 3
137graph_color <- "darkseagreen4"
138geom_thickness <- 0.25
139geom_age_orientation <- coord_flip()
140ENABLE_LOG_SCALE <- TRUE
141
142### If TRUE, graph total cases, new cases, and total deaths on a log scale
143### If FALSE, graph them normally
144if (isTRUE(ENABLE_LOG_SCALE)) {
145 graph_log10_opts1 <- scale_y_log10(breaks = c(1, 10, 100, 1000))
146 graph_log10_opts2 <- annotation_logticks()
147} else {
148 graph_log10_opts1 <- geom_blank()
149 graph_log10_opts2 <- geom_blank()
150}
151
152
153################################################################################
154### Adjusting the geometry of our map
155################################################################################
156### Massage the state/county map data a bit
157state_mapdata <- tn_pop_df
158
159### Make a copy of the state geometry
160counties <- state_mapdata
161
162### Find the center of each county so we can add the number of infected
163county_centers <-
164 tn_map_df$geometry %>%
165 as_Spatial() %>%
166 gCentroid(byid = TRUE)
167
168### The centroid function is very good, but occasionally can give less-than-
169### perfect results, often due to counties with extreme concave shape. We
170### can use "nudges" to adjust the location of the number to make it look good.
171### Default nudge is 0.
172counties$nudge_y <- 0
173counties$nudge_x <- 0
174
175counties$nudge_y[counties$County == "Chester"] <- 0.035
176counties$nudge_y[counties$County == "Houston"] <- 0.015
177counties$nudge_y[counties$County == "Putnam"] <- 0.025
178counties$nudge_y[counties$County == "Loudon"] <- 0.015
179counties$nudge_y[counties$County == "Anderson"] <- -0.015
180counties$nudge_x[counties$County == "Pickett"] <- -0.065
181counties$nudge_x[counties$County == "Unicoi"] <- -0.070
182counties$nudge_y[counties$County == "Unicoi"] <- -0.035
183counties$nudge_x[counties$County == "Hancock"] <- -0.035
184counties$nudge_x[counties$County == "Trousdale"] <- -0.020
185counties$nudge_x[counties$County == "Lake"] <- 0.030
186counties$nudge_x[counties$County == "Moore"] <- 0.005
187counties$nudge_y[counties$County == "Moore"] <- 0.015
188
189
190################################################################################
191### Math section for maps
192################################################################################
193### Get the number of total confirmed cases by county
194totalinfected_by_county <-
195 infected_df %>%
196 arrange(Date) %>%
197 tail(n = 1) %>%
198 select(-all_of(non_county_fields)) %>%
199 gather(Date) %>%
200 rename(County = Date) %>%
201 rename(totalinfected = value)
202
203### Get the number of new confirmed cases by county
204new_date <-
205 infected_df %>%
206 pull(Date) %>%
207 max() %>%
208 as_tibble() %>%
209 rename(Date = value) %>%
210 pull()
211
212newinfected_by_county <-
213 infected_df %>%
214 arrange(Date) %>%
215 tail(n = 2) %>%
216 select(-all_of(non_county_fields))
217
218newinfected_by_county <-
219 (newinfected_by_county[2, ] - newinfected_by_county[1, ]) %>%
220 as_tibble() %>%
221 add_column(new_date, .before = 1) %>%
222 rename(Date = new_date) %>%
223 gather(Date) %>%
224 rename(County = Date) %>%
225 rename(newinfected = value)
226
227### Get the number of total deaths by county
228totaldeaths_by_county <-
229 deaths_df %>%
230 filter(Date > "2020-03-31") %>%
231 head(n = 1) %>%
232 select("Date", all_of(all_counties)) %>%
233 gather(Date) %>%
234 rename(County = Date) %>%
235 rename(totaldeaths = value)
236
237### Get the number of new deaths by county
238new_date <-
239 deaths_df %>%
240 select("Date") %>%
241 filter(!is.na(Date)) %>%
242 pull(Date) %>%
243 max() %>%
244 as_tibble() %>%
245 pull()
246
247newdeaths_by_county <-
248 deaths_df %>%
249 filter(!is.na(Date)) %>%
250 arrange(Date) %>%
251 tail(n = 2) %>%
252 select(-all_of(non_county_fields))
253
254newdeaths_by_county <-
255 (newdeaths_by_county[2, ] - newdeaths_by_county[1, ]) %>%
256 as_tibble() %>%
257 add_column(new_date, .before = 1) %>%
258 rename(Date = new_date) %>%
259 gather(Date) %>%
260 rename(County = Date) %>%
261 rename(newdeaths = value)
262
263### Get the number of total deaths by county
264recovered_by_county <-
265 recover_df %>%
266 filter(Date > "2020-03-31") %>%
267 arrange(Date) %>%
268 tail(n = 1) %>%
269 select("Date", all_of(all_counties)) %>%
270 gather(Date) %>%
271 rename(County = Date) %>%
272 rename(recovered = value)
273
274### Still sick
275sick_by_county <-
276 totalinfected_by_county %>%
277 left_join(recovered_by_county, by = "County") %>%
278 left_join(totaldeaths_by_county, by = "County") %>%
279 mutate(sick = totalinfected - recovered - totaldeaths) %>%
280 select(County, sick)
281
282### Combine the data
283stats_df <-
284 totalinfected_by_county %>%
285 left_join(newinfected_by_county, by = "County") %>%
286 left_join(totaldeaths_by_county, by = "County") %>%
287 left_join(newdeaths_by_county, by = "County") %>%
288 left_join(sick_by_county, by = "County") %>%
289 left_join(recovered_by_county, by = "County")
290
291### Combine all the above together
292counties <-
293 counties %>%
294 left_join(tn_map_df, by = "GEOID") %>%
295 left_join(stats_df, by = "County") %>%
296 mutate(newinfected = if_else(newinfected == -1, 0, newinfected))
297
298
299################################################################################
300### Map of total COVID-19 cases
301################################################################################
302### We want to add text with the number of infected in each county to that
303### county on the map. But we *don't* want to label counties that don't yet
304### have any confirmed cases.
305totalinfected_label <- counties$totalinfected
306totalinfected_label[totalinfected_label == 0] <- ""
307
308### By default use black for the font, but if the value is over 1/2's of the way
309### to the max, use the white font instead as it looks nicer
310frac <- 0.5 * max(counties$totalinfected)
311counties$textcolor <- if_else(counties$totalinfected >= frac, "white", "black")
312
313### Map: COVID-19 confirmed infected per county
314p_map_totalinfected <-
315 ggplot(counties) +
316 theme_void() +
317 theme(legend.title = element_blank()) +
318 theme(plot.title = element_text(hjust = 0.5)) +
319 theme(legend.position = "none") +
320 geom_sf(data = counties$geometry, aes(fill = counties$totalinfected),
321 size = geom_thickness) +
322 geom_text(data = counties, size = map_fontsize, fontface = "bold",
323 color = counties$textcolor,
324 aes(x = county_centers$x,
325 y = county_centers$y,
326 label = totalinfected_label),
327 nudge_x = counties$nudge_x,
328 nudge_y = counties$nudge_y) +
329 labs(title = "Total Cases") +
330 scale_fill_gradientn(colours = map_palette)
331
332
333################################################################################
334### Map of new COVID-19 cases
335################################################################################
336### We want to add text with the number of infected in each county to that
337### county on the map. But we *don't* want to label counties that don't yet
338### have any confirmed cases.
339newinfected_label <- counties$newinfected
340newinfected_label[newinfected_label <= 0] <- ""
341
342### By default use black for the font, but if the value is over 1/2's of the way
343### to the max, use the white font instead as it looks nicer
344frac <- 0.5 * max(counties$newinfected)
345counties$textcolor <- if_else(counties$newinfected >= frac, "white", "black")
346
347### Map: COVID-19 confirmed new infected per county
348p_map_newinfected <-
349 ggplot(counties) +
350 theme_void() +
351 theme(legend.title = element_blank()) +
352 theme(plot.title = element_text(hjust = 0.5)) +
353 theme(legend.position = "none") +
354 geom_sf(data = counties$geometry, aes(fill = counties$newinfected),
355 size = geom_thickness) +
356 geom_text(data = counties, size = map_fontsize, fontface = "bold",
357 color = counties$textcolor,
358 aes(x = county_centers$x,
359 y = county_centers$y,
360 label = newinfected_label),
361 nudge_x = counties$nudge_x,
362 nudge_y = counties$nudge_y) +
363 labs(title = "New Cases") +
364 scale_fill_gradientn(colours = map_palette)
365
366
367################################################################################
368### Map of deaths by county
369################################################################################
370total_deaths_label <- counties$totaldeaths
371total_deaths_label[total_deaths_label == 0] <- ""
372
373### By default use black for the font, but if the value is over 1/2's of the way
374### to the max, use the white font instead as it looks nicer
375frac <- 0.5 * max(counties$totaldeaths)
376counties$textcolor <- if_else(counties$totaldeaths >= frac, "white", "black")
377
378p_map_totaldeaths <- ggplot(counties) +
379 theme_void() +
380 theme(legend.title = element_blank()) +
381 theme(plot.title = element_text(hjust = 0.5)) +
382 theme(legend.position = "none") +
383 geom_sf(data = counties$geometry, size = geom_thickness,
384 aes(fill = counties$totaldeaths)) +
385 geom_text(data = counties, size = map_fontsize, fontface = "bold",
386 color = counties$textcolor,
387 aes(x = county_centers$x,
388 y = county_centers$y,
389 label = total_deaths_label),
390 nudge_x = counties$nudge_x,
391 nudge_y = counties$nudge_y) +
392 labs(title = "Total Deaths") +
393 scale_fill_gradientn(colours = map_palette)
394
395new_deaths_label <- counties$newdeaths
396new_deaths_label[new_deaths_label == 0] <- ""
397
398### By default use black for the font, but if the value is over 1/2's of the way
399### to the max, use the white font instead as it looks nicer
400frac <- 0.5 * max(counties$newdeaths)
401counties$textcolor <- if_else(counties$newdeaths > frac, "white", "black")
402
403p_map_newdeaths <- ggplot(counties) +
404 theme_void() +
405 theme(legend.title = element_blank()) +
406 theme(plot.title = element_text(hjust = 0.5)) +
407 theme(legend.position = "none") +
408 geom_sf(data = counties$geometry, aes(fill = counties$newdeaths),
409 size = geom_thickness) +
410 geom_text(data = counties, size = map_fontsize, fontface = "bold",
411 color = counties$textcolor,
412 aes(x = county_centers$x,
413 y = county_centers$y,
414 label = new_deaths_label),
415 nudge_x = counties$nudge_x,
416 nudge_y = counties$nudge_y) +
417 labs(title = "New Deaths") +
418 scale_fill_gradientn(colours = map_palette)
419
420
421################################################################################
422### Map of deaths per 100k by county
423################################################################################
424counties$deaths_per100k <- round(100000 * counties$totaldeaths /
425 counties$POP2018, 3)
426deaths100k_label <- round(counties$deaths_per100k, 1)
427deaths100k_label[deaths100k_label == 0] <- ""
428
429### By default use black for the font, but if the value is over 1/2's of the way
430### to the max, use the white font instead as it looks nicer
431frac <- 0.5 * max(counties$deaths_per100k)
432counties$textcolor <- if_else(counties$deaths_per100k >= frac, "white", "black")
433
434p_map_totaldeaths_100k <- ggplot(counties) +
435 theme_void() +
436 theme(legend.title = element_blank()) +
437 theme(plot.title = element_text(hjust = 0.5)) +
438 theme(legend.position = "none") +
439 geom_sf(data = counties$geometry, aes(fill = counties$deaths_per100k),
440 size = geom_thickness) +
441 geom_text(data = counties, size = map_fontsize, fontface = "bold",
442 color = counties$textcolor,
443 aes(x = county_centers$x,
444 y = county_centers$y,
445 label = deaths100k_label),
446 nudge_x = counties$nudge_x,
447 nudge_y = counties$nudge_y) +
448 labs(title = "Total Deaths per 100k") +
449 scale_fill_gradientn(colours = map_palette)
450
451
452################################################################################
453### Map of total COVID-19 cases per 100k
454################################################################################
455counties$infected_percapita <-
456 round(100000 * counties$totalinfected / counties$POP2018, 0)
457
458frac <- 0.5 * max(counties$infected_percapita)
459counties$textcolor <- if_else(counties$infected_percapita >= frac,
460 "white", "black")
461
462infected_label_per <- round(counties$infected_percapita, 1)
463infected_label_per[infected_label_per == 0] <- ""
464
465p_map_totalinfected_per <- ggplot(counties) +
466 theme_void() +
467 theme(plot.title = element_text(hjust = 0.5)) +
468 theme(plot.caption = element_text(hjust = 0.5)) +
469 theme(legend.title = element_blank()) +
470 theme(legend.position = "none") +
471 geom_sf(data = counties$geometry, size = geom_thickness,
472 aes(fill = counties$infected_percapita)) +
473 geom_text(data = counties, size = map_fontsize, fontface = "bold",
474 color = counties$textcolor,
475 aes(x = county_centers$x,
476 y = county_centers$y,
477 label = infected_label_per),
478 nudge_x = counties$nudge_x,
479 nudge_y = counties$nudge_y) +
480 labs(title = "Total Cases per 100k") +
481 scale_fill_gradientn(colours = map_palette)
482
483
484################################################################################
485### Map of people still sick (not recovered) by county
486################################################################################
487frac <- 0.5 * max(counties$sick)
488counties$textcolor <- if_else(counties$sick >= frac, "white", "black")
489
490sick_label_per <- counties$sick
491sick_label_per[sick_label_per == 0] <- ""
492
493p_map_sick <- ggplot(counties) +
494 theme_void() +
495 theme(plot.title = element_text(hjust = 0.5)) +
496 theme(plot.caption = element_text(hjust = 0.5)) +
497 theme(legend.title = element_blank()) +
498 theme(legend.position = "none") +
499 geom_sf(data = counties$geometry, size = geom_thickness,
500 aes(fill = counties$sick)) +
501 geom_text(data = counties, size = map_fontsize, fontface = "bold",
502 color = counties$textcolor,
503 aes(x = county_centers$x,
504 y = county_centers$y,
505 label = sick_label_per),
506 nudge_x = counties$nudge_x,
507 nudge_y = counties$nudge_y) +
508 labs(title = "Still Sick") +
509 scale_fill_gradientn(colours = map_palette)
510
511
512################################################################################
513### Map of people recovered by county
514################################################################################
515frac <- 0.5 * max(counties$recovered)
516counties$textcolor <- if_else(counties$recovered >= frac, "white", "black")
517
518recovered_label_per <- counties$recovered
519recovered_label_per[recovered_label_per == 0] <- ""
520
521p_map_recovered <- ggplot(counties) +
522 theme_void() +
523 theme(plot.title = element_text(hjust = 0.5)) +
524 theme(plot.caption = element_text(hjust = 0.5)) +
525 theme(legend.title = element_blank()) +
526 theme(legend.position = "none") +
527 geom_sf(data = counties$geometry, size = geom_thickness,
528 aes(fill = counties$recovered)) +
529 geom_text(data = counties, size = map_fontsize, fontface = "bold",
530 color = counties$textcolor,
531 aes(x = county_centers$x,
532 y = county_centers$y,
533 label = recovered_label_per),
534 nudge_x = counties$nudge_x,
535 nudge_y = counties$nudge_y) +
536 labs(title = "Recovered") +
537 scale_fill_gradientn(colours = map_palette)
538
539
540################################################################################
541### Math section for charts
542################################################################################
543infected_data <-
544 infected_df %>%
545 arrange(Date) %>%
546 select("Date", all_of(graph_counties)) %>%
547 mutate(Total = rowSums(select(., !starts_with("Date")))) %>%
548 select(Date, Total) %>%
549 rename(totalinfected = Total) %>%
550 rename(date = Date) %>%
551 filter(totalinfected != 0) %>%
552 mutate(newinfected = c(totalinfected[1], diff(totalinfected))) %>%
553 mutate(numdays = difftime(as.POSIXct(date),
554 as.POSIXct(as.Date(min(date)) - 1, tz = "UTC"),
555 units = "days")) %>%
556 mutate(numdays = as.numeric(numdays))
557
558hospital_data <-
559 hospital_df %>%
560 arrange(Date) %>%
561 mutate(numdays = difftime(as.POSIXct(Date),
562 as.POSIXct(as.Date(min(Date)) - 1, tz = "UTC"),
563 units = "days")) %>%
564 mutate(numdays = as.numeric(numdays)) %>%
565 rename(total_hospitalizations = Hospitalizations) %>%
566 mutate(new_hospitalizations = c(total_hospitalizations[1],
567 diff(total_hospitalizations)))
568
569recovered_data <-
570 recover_df %>%
571 rename(Total = Recovered) %>%
572 select("Date", all_of(graph_counties)) %>%
573 mutate(Total = rowSums(select(., !starts_with("Date")))) %>%
574 select(Date, Total) %>%
575 rename(totalrecovered = Total) %>%
576 rename(date = Date) %>%
577 filter(totalrecovered != 0) %>%
578 arrange(date) %>%
579 mutate(newrecovered = c(totalrecovered[1], diff(totalrecovered))) %>%
580 mutate(numdays = difftime(as.POSIXct(date),
581 as.POSIXct(as.Date(min(date)) - 1, tz = "UTC"),
582 units = "days")) %>%
583 mutate(numdays = as.numeric(numdays))
584
585deaths_data <-
586 deaths_df %>%
587 na.omit() %>%
588 arrange(Date) %>%
589 select("Date", all_of(graph_counties)) %>%
590 mutate(Total = rowSums(select(., !starts_with("Date")))) %>%
591 mutate(New = c(1, diff(Total))) %>%
592 select(Date, Total, New) %>%
593 filter(Total != 0) %>%
594 rename(totaldeaths = Total) %>%
595 rename(newdeaths = New) %>%
596 rename(date = Date) %>%
597 mutate(numdays = difftime(as.POSIXct(date),
598 as.POSIXct(as.Date(min(date)) - 1, tz = "UTC"),
599 units = "days")) %>%
600 mutate(numdays = as.numeric(numdays))
601
602sick_data <-
603 infected_data %>%
604 left_join(recovered_data, by = "date") %>%
605 left_join(deaths_data, by = "date") %>%
606 select(date, totalinfected, totaldeaths, totalrecovered) %>%
607 filter(date >= as.Date("2020-04-10")) %>%
608 filter(!is.na(totaldeaths)) %>%
609 filter(!is.na(totalrecovered)) %>%
610 filter(!is.na(totalinfected)) %>%
611 mutate(current_sick = totalinfected - totaldeaths - totalrecovered) %>%
612 select(date, current_sick)
613
614
615################################################################################
616### Graph: COVID-19 Total Confirmed Cases
617################################################################################
618### Use non-linear fitting to do the exponential fit
619fit_totalinfected <-
620 nls(totalinfected ~ a * exp(b * numdays),
621 data = infected_data,
622 start = list(a = 334,
623 b = 0.0677))
624
625### Extract the regression coefficents so we can add the regression formula
626### to the graph
627intercept <- coef(fit_totalinfected)[1]
628d_constant <- coef(fit_totalinfected)[2]
629
630### Add the fit to our dataframe
631infected_data <-
632 infected_data %>%
633 mutate(fit_totalinfected = intercept *
634 exp(d_constant * infected_data$numdays))
635
636### Compute the doubling time
637doubling_time <- log(2) / d_constant
638
639### Compute the "order-of-magnitude" time
640magnitude_time <- log(10) / d_constant
641
642### Growth rate
643growth_rate <- exp(log(2) / doubling_time) - 1
644
645### What's the current number of total infected cases
646total_num <-
647 infected_data %>%
648 tail(n = 1) %>%
649 select("totalinfected") %>%
650 pull()
651
652### What's the date of the last data point
653data_date <- infected_data %>% select("date") %>% tail(n = 1) %>% pull()
654
655p_total <- ggplot(data = infected_data) +
656 theme_classic() +
657 theme(axis.text.x = element_text(angle = 45)) +
658 theme(axis.text.x = element_text(vjust = 0.7)) +
659 theme(axis.text.x = element_text(hjust = 0.8)) +
660 theme(legend.title = element_blank()) +
661 theme(legend.position = "none") +
662 geom_point(data = infected_data, shape = 19, size = 0.5,
663 aes(x = as.Date(date), y = totalinfected), color = "black") +
664 geom_line(data = infected_data, color = graph_color,
665 aes(x = as.Date(date), y = fit_totalinfected)) +
666 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
667 graph_log10_opts1 +
668 graph_log10_opts2 +
669 labs(title = paste("Total Cases: ", total_num, sep = ""), x = "", y = "")
670
671
672################################################################################
673### Graph: COVID-19 New Confirmed Cases
674################################################################################
675### If we assume infected = a * B^(c * days), then
676###
677### d(infected)/dt = c * log_e(B) * infected
678###
679### We can use the same constants from the total fit on the new fit. If new
680### cases are still growing exponentially, it will work. And when it doesn't
681### work, that's a sign that we're either "flattening the curve", or something
682### else weird is going on
683
684### Add the fit to our dataframe
685infected_data <-
686 infected_data %>%
687 mutate(fit_newinfected = fit_totalinfected * d_constant)
688
689### What's the current number of total infected cases
690new_num <- infected_data %>% tail(n = 1) %>% select("newinfected") %>% pull()
691SMA_num <- infected_data$newinfected %>% SMA(n = 7) %>% tail(n = 1) %>% round()
692
693### Graph: new COVID-19 Cases
694p_new <- ggplot(data = infected_data) +
695 theme_classic() +
696 theme(axis.text.x = element_text(angle = 45)) +
697 theme(axis.text.x = element_text(vjust = 0.7)) +
698 theme(axis.text.x = element_text(hjust = 0.8)) +
699 theme(legend.title = element_blank()) +
700 theme(legend.position = "none") +
701 geom_point(data = infected_data, shape = 19, size = 0.5,
702 aes(x = as.Date(date), y = newinfected), color = "black") +
703 geom_line(data = infected_data, color = graph_color,
704 size = 1.0, aes(x = as.Date(date) - 3,
705 y = SMA(newinfected, n = 7))) +
706 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
707 graph_log10_opts1 +
708 graph_log10_opts2 +
709 labs(title = paste("New Cases: ", new_num, ", SMA: ", SMA_num, sep = ""),
710 x = "", y = "")
711
712
713################################################################################
714### Graph: COVID-19 Cases [ log(Total) vs log(New) ]
715### Reference: https://aatishb.com/covidtrends/
716################################################################################
717p_total_new <- ggplot(data = infected_data,
718 aes(x = totalinfected, y = newinfected)) +
719 theme_classic() +
720 theme(axis.text.x = element_text(angle = 45)) +
721 theme(axis.text.x = element_text(vjust = 0.7)) +
722 theme(axis.text.x = element_text(hjust = 0.8)) +
723 theme(legend.position = "none") +
724 geom_point(data = infected_data, shape = 19, size = 0.5,
725 aes(x = totalinfected, y = newinfected), color = "black") +
726 geom_smooth(method = "loess", formula = y ~ x,
727 se = FALSE, size = 0.5, color = graph_color) +
728 scale_x_log10(breaks = c(1, 10, 100, 1000)) +
729 scale_y_log10(breaks = c(1, 10, 100, 1000)) +
730 annotation_logticks() +
731 labs(title = "Log(New vs Total)", x = "", y = "")
732
733
734################################################################################
735### Graph: COVID-19 Hospitalizations
736################################################################################
737hospital_num <-
738 hospital_data %>%
739 arrange(Date) %>%
740 tail(n = 1) %>%
741 select("total_hospitalizations") %>%
742 pull()
743
744p_total_hospital <- ggplot(data = hospital_data,
745 aes(x = as.Date(Date), y = total_hospitalizations)) +
746 theme_classic() +
747 theme(axis.text.x = element_text(angle = 45)) +
748 theme(axis.text.x = element_text(vjust = 0.7)) +
749 theme(axis.text.x = element_text(hjust = 0.8)) +
750 theme(legend.title = element_blank()) +
751 theme(legend.position = "none") +
752 geom_bar(stat = "identity", fill = graph_color) +
753 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
754 labs(title = paste("TN Hospitalizations: ", hospital_num, sep = ""),
755 x = "", y = "")
756
757
758################################################################################
759### Graph: COVID-19 Total Recovered
760################################################################################\
761new_num <-
762 recovered_data %>%
763 tail(n = 1) %>%
764 select("totalrecovered") %>%
765 pull()
766
767p_recover <- ggplot(data = recovered_data,
768 aes(x = as.Date(date), y = totalrecovered)) +
769 theme_classic() +
770 theme(axis.text.x = element_text(angle = 45)) +
771 theme(axis.text.x = element_text(vjust = 0.7)) +
772 theme(axis.text.x = element_text(hjust = 0.8)) +
773 theme(legend.title = element_blank()) +
774 theme(legend.position = "none") +
775 geom_point(data = recovered_data, shape = 19, size = 0.5,
776 aes(x = as.Date(date), y = totalrecovered), color = "black") +
777 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
778 labs(title = paste("Total Recovered: ", new_num, sep = ""), x = "", y = "")
779
780
781################################################################################
782### Graph: COVID-19 Total Infected by Age
783################################################################################
784p_age_infected <- ggplot(data = age_df, aes(x = infected_statewide, y = age)) +
785 theme_classic() +
786 theme(axis.text.x = element_text(angle = 45)) +
787 theme(axis.text.x = element_text(vjust = 0.7)) +
788 theme(axis.text.x = element_text(hjust = 0.8)) +
789 geom_age_orientation +
790 geom_bar(stat = "identity", fill = graph_color) +
791 labs(title = "TN Infected by Age", x = "", y = "")
792
793
794################################################################################
795### Graph: COVID-19 Total Infected by Age
796################################################################################
797p_age_deaths <- ggplot(data = age_df, aes(x = deaths_statewide, y = age)) +
798 theme_classic() +
799 theme(axis.text.x = element_text(angle = 45)) +
800 theme(axis.text.x = element_text(vjust = 0.7)) +
801 theme(axis.text.x = element_text(hjust = 0.8)) +
802 geom_age_orientation +
803 geom_bar(stat = "identity", fill = graph_color) +
804 labs(title = "TN Deaths by Age", x = "", y = "")
805
806
807################################################################################
808### Graph: COVID-19 Total/New Deaths
809################################################################################
810### Use non-linear fitting to do the exponential fit
811fit_totaldeaths <- nls(totaldeaths ~ a * exp(b * numdays),
812 data = deaths_data, start = list(a = 12.51, b = 0.082))
813
814### Extract the regression coefficents so we can add the regression formula
815### to the graph
816intercept <- coef(fit_totaldeaths)[1]
817d_constant <- coef(fit_totaldeaths)[2]
818
819### Add fit for total/new deaths
820deaths_data <-
821 deaths_data %>%
822 mutate(fit_totaldeaths = intercept * exp(d_constant * numdays)) %>%
823 mutate(fit_newdeaths = fit_totaldeaths * d_constant)
824
825p_totaldeaths <- ggplot(data = deaths_data) +
826 theme_classic() +
827 theme(axis.text.x = element_text(angle = 45)) +
828 theme(axis.text.x = element_text(vjust = 0.7)) +
829 theme(axis.text.x = element_text(hjust = 0.8)) +
830 theme(legend.title = element_blank()) +
831 theme(legend.position = "none") +
832 geom_point(shape = 19, size = 0.5,
833 aes(x = as.Date(date), y = totaldeaths), color = "black") +
834 geom_line(color = graph_color,
835 aes(x = as.Date(date), y = fit_totaldeaths)) +
836 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
837 graph_log10_opts1 +
838 graph_log10_opts2 +
839 labs(title = paste("Total Deaths: ", max(deaths_data$totaldeaths), sep = ""),
840 x = "", y = "")
841
842p_newdeaths <- ggplot(data = deaths_data) +
843 theme_classic() +
844 theme(axis.text.x = element_text(angle = 45)) +
845 theme(axis.text.x = element_text(vjust = 0.7)) +
846 theme(axis.text.x = element_text(hjust = 0.8)) +
847 theme(legend.title = element_blank()) +
848 theme(legend.position = "none") +
849 geom_point(shape = 19, size = 0.5,
850 aes(x = as.Date(date), y = newdeaths), color = "black") +
851 geom_line(color = graph_color,
852 aes(x = as.Date(date), y = fit_newdeaths)) +
853 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
854 graph_log10_opts1 +
855 graph_log10_opts2 +
856 labs(title = paste("New Deaths: ", max(deaths_data$newdeaths), sep = ""),
857 x = "", y = "")
858
859
860################################################################################
861### Graph: COVID-19 Current Sick
862################################################################################
863new_num <- sick_data %>% tail(n = 1) %>% select("current_sick") %>% pull()
864
865p_sick <- ggplot(data = sick_data, aes(x = as.Date(date), y = current_sick)) +
866 theme_classic() +
867 theme(axis.text.x = element_text(angle = 45)) +
868 theme(axis.text.x = element_text(vjust = 0.7)) +
869 theme(axis.text.x = element_text(hjust = 0.8)) +
870 theme(legend.title = element_blank()) +
871 theme(legend.position = "none") +
872 geom_point(data = sick_data, shape = 19, size = 0.5,
873 aes(x = as.Date(date), y = current_sick), color = "black") +
874 scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
875 labs(title = paste("Currently Sick: ", new_num, sep = ""), x = "", y = "")
876
877
878################################################################################
879### Arrange the graphs into the final product
880################################################################################
881# What's the date of the last data point
882data_date <- infected_df %>% arrange(Date) %>% tail(n = 1) %>% pull("Date")
883
884# Title
885title_string <- paste("COVID-19 Confirmed Cases in ", location,
886 " [", as.Date(data_date), "]", sep = "")
887
888# Footer
889footer_string <- paste("Source: Tennessee Department of Health\n",
890 "https://www.tn.gov/health/cedep/ncov.html", sep = "")
891
892grob_charts <-
893 plot_grid(p_total, p_new, p_totaldeaths, p_recover, p_sick,
894 p_age_infected, p_age_deaths, p_total_hospital, p_total_new, NA,
895 ncol = 5, nrow = 2, align = "hv")
896
897grob_maps <-
898 plot_grid(p_map_recovered, p_map_totalinfected,
899 p_map_sick, p_map_totalinfected_per,
900 p_map_totaldeaths, p_map_newinfected,
901 ncol = 2, nrow = 3, align = "hv")
902
903charts <- plot_grid(grob_maps,
904 grob_charts,
905 ncol = 1, nrow = 2,
906 rel_widths = c(1.2, 1),
907 rel_heights = c(1, 1))
908
909title <- ggdraw() + draw_label(title_string, fontface = "bold")
910Sys.sleep(2)
911
912footer <- ggdraw() + draw_label(footer_string, size = 10)
913Sys.sleep(2)
914
915final <- plot_grid(title, charts, footer, nrow = 3,
916 rel_heights = c(0.05, 1, 0.05))
917Sys.sleep(2)
918
919print(final)
920
921### Save image to file and we're done!
922scale <- 1.3
923width <- round(1300 * scale)
924height <- round(900 * scale)
925png("TN_COVID_Viridis.png", width = width, height = height)
926plot(final)
927dev.off()