TN COVID-19 script, Apr 21

################################################################################
### Render a map of Tennessee using TN Dept. of Health data on COVID-19
### infections showing total cases, new cases, cases by age, and a map of
### confirmed infected.
###
### TN Dept. of Health COVID-19 website:
###     https://www.tn.gov/health/cedep/ncov.html
###
### Note: You will need a Census API key to access their data.  You
###       can get a key by going to this website:
###
###       https://api.census.gov/data/key_signup.html
###
### Once you have the key, just add it below.
################################################################################

### Set the Census API key
census_key <- "PUT_YOUR_CENSUS_API_KEY_HERE"

### Let's start by cleaning the environment, it currently causes the graphs
### some problems if the environment is already populated
rm(list = ls())

# Need to install the github version of googlesheets4
# install.packages("devtools")
devtools::install_github("tidyverse/googlesheets4")

### Load the necessary libraries
packages <- c("tidyverse", "dplyr", "readxl", "ggplot2", "sf", "rgeos", "TTR",
              "scales", "cowplot", "viridis", "gridExtra", "tidycensus",
              "googlesheets4")

new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages)) install.packages(new_packages, quiet = TRUE)
invisible(lapply(packages, "library", quietly = TRUE,
                 character.only = TRUE, warn.conflicts = FALSE))

### Set the Census API key
census_api_key(census_key)

### Enable map caching
options(tigris_use_cache = TRUE)


################################################################################
### Read data from my spreadsheet
################################################################################
### Load COVID data from Google Sheets
gs4_deauth()
covid_data <- "https://docs.google.com/spreadsheets/d/1Opqs3DIYTlPbPhqrUa66RfKBBDVo0u5sntFM6uYZDzk/edit?usp=sharing"
infected_df  <- read_sheet(covid_data, sheet = "Infected")
deaths_df    <- read_sheet(covid_data, sheet = "Deaths")
hospital_df  <- read_sheet(covid_data, sheet = "Hospitalizations")
testres_df   <- read_sheet(covid_data, sheet = "Test_Results")
testrates_df <- read_sheet(covid_data, sheet = "Testing_Rates") %>%
                mutate(County = gsub(" County", "", County))
recover_df   <- read_sheet(covid_data, sheet = "Recovered")
age_df       <- read_sheet(covid_data, sheet = "Age_Data")


################################################################################
### Subset the the data here (if desired)
################################################################################
### For the map, we only want the latest total for each county.
non_county_fields <- c("Date", "Total", "Other_State_Country",
                       "Unknown", "ERROR!!!")

### All counties
all_counties <-
  colnames(infected_df) %>%
  enframe() %>%
  filter(!colnames(infected_df) %in% non_county_fields) %>%
  select("value") %>%
  deframe()

### Counties in the Nashville Metropolitan Statistical Area
nashville_msa <-
  c("Cannon", "Cheatham", "Davidson", "Dickson", "Macon", "Maury", "Robertson",
    "Rutherford", "Smith", "Sumner", "Trousdale", "Williamson", "Wilson")

middle_tn <- c(
  "Bedford",    "Cannon",  "Cheatham",  "Clay",       "Coffee",    "Davidson",
  "DeKalb",     "Dickson", "Fentress",  "Franklin",   "Giles",     "Grundy",
  "Hickman",    "Houston", "Humphreys", "Jackson",    "Lawrence",  "Lewis",
  "Lincoln",    "Macon",   "Marshall",  "Maury",      "Montgomery", "Moore",
  "Overton",    "Perry",   "Pickett",   "Putnam",     "Robertson", "Rutherford",
  "Sequatchie", "Smith",   "Stewart",   "Sumner",     "Trousdale", "Van Buren",
  "Warren",     "Wayne",   "White",     "Williamson", "Wilson")


################################################################################
### Select the map you want
################################################################################
#graph_counties <- nashville_msa
#map_counties <- nashville_msa
#location <- "Nashville MSA"

#graph_counties <- c("Davidson")
#map_counties <- c("Davidson")
#location <- "Davidson County"

#graph_counties <- middle_tn
#map_counties <- middle_tn
#location <- "Middle TN"

graph_counties <- c("Total")
map_counties <- all_counties
location <- "Tennessee"

### Use TidyCensus to pull population data as well as map geometry
### from the Census bureau
county_acs <-
  get_acs(geography   = "county",
          variables   = c(POP2018 = "B01003_001"),
          state       = c("47"),
          county      = map_counties,
          year        = 2018,
          geometry    = TRUE,
          cache_table = TRUE) %>%
  rename(POP2018 = estimate)

### Split the Census data apart into a map...
tn_map_df <-
  county_acs %>%
  select(GEOID, geometry)

### And a population dataframe
tn_pop_df <-
  county_acs %>%
  st_set_geometry(NULL) %>%
  mutate(NAME = str_replace(NAME, " County, Tennessee", "")) %>%
  select(GEOID, NAME, POP2018) %>%
  rename(County = NAME)


################################################################################
### A few visual configuration settings
################################################################################
map_palette <- c("#c4c4c4", viridis(n = 1000, direction = -1))
map_fontsize <- 3
graph_color <- "darkseagreen4"
geom_thickness <- 0.25
geom_age_orientation <- coord_flip()
ENABLE_LOG_SCALE <- TRUE

### If TRUE, graph total cases, new cases, and total deaths on a log scale
### If FALSE, graph them normally
if (isTRUE(ENABLE_LOG_SCALE)) {
     graph_log10_opts1 <- scale_y_log10(breaks = c(1, 10, 100, 1000))
     graph_log10_opts2 <- annotation_logticks()
} else {
     graph_log10_opts1 <- geom_blank()
     graph_log10_opts2 <- geom_blank()
}


################################################################################
### Adjusting the geometry of our map
################################################################################
### Massage the state/county map data a bit
state_mapdata <- tn_pop_df

### Make a copy of the state geometry
counties <- state_mapdata

### Find the center of each county so we can add the number of infected
county_centers <-
  tn_map_df$geometry %>%
  as_Spatial() %>%
  gCentroid(byid = TRUE)

### The centroid function is very good, but occasionally can give less-than-
### perfect results, often due to counties with extreme concave shape.   We
### can use "nudges" to adjust the location of the number to make it look good.
### Default nudge is 0.
counties$nudge_y <- 0
counties$nudge_x <- 0

counties$nudge_y[counties$County == "Chester"]   <-  0.035
counties$nudge_y[counties$County == "Houston"]   <-  0.015
counties$nudge_y[counties$County == "Putnam"]    <-  0.025
counties$nudge_y[counties$County == "Loudon"]    <-  0.015
counties$nudge_y[counties$County == "Anderson"]  <- -0.015
counties$nudge_x[counties$County == "Pickett"]   <- -0.065
counties$nudge_x[counties$County == "Unicoi"]    <- -0.070
counties$nudge_y[counties$County == "Unicoi"]    <- -0.035
counties$nudge_x[counties$County == "Hancock"]   <- -0.035
counties$nudge_x[counties$County == "Trousdale"] <- -0.020
counties$nudge_x[counties$County == "Lake"]      <-  0.030
counties$nudge_x[counties$County == "Moore"]     <-  0.005
counties$nudge_y[counties$County == "Moore"]     <-  0.015


################################################################################
### Math section for maps
################################################################################
### Get the number of total confirmed cases by county
totalinfected_by_county <-
  infected_df %>%
  arrange(Date) %>%
  tail(n = 1) %>%
  select(-all_of(non_county_fields)) %>%
  gather(Date) %>%
  rename(County = Date) %>%
  rename(totalinfected = value)

### Get the number of new confirmed cases by county
new_date <-
  infected_df %>%
  pull(Date) %>%
  max() %>%
  as_tibble() %>%
  rename(Date = value) %>%
  pull()

newinfected_by_county <-
  infected_df %>%
  arrange(Date) %>%
  tail(n = 2) %>%
  select(-all_of(non_county_fields))

newinfected_by_county <-
  (newinfected_by_county[2, ] - newinfected_by_county[1, ]) %>%
  as_tibble() %>%
  add_column(new_date, .before = 1) %>%
  rename(Date = new_date) %>%
  gather(Date) %>%
  rename(County = Date) %>%
  rename(newinfected = value)

### Get the number of total deaths by county
totaldeaths_by_county <-
  deaths_df %>%
  filter(Date > "2020-03-31") %>%
  head(n = 1) %>%
  select("Date", all_of(all_counties)) %>%
  gather(Date) %>%
  rename(County = Date) %>%
  rename(totaldeaths = value)

### Get the number of new deaths by county
new_date <-
  deaths_df %>%
  select("Date") %>%
  filter(!is.na(Date)) %>%
  pull(Date) %>%
  max() %>%
  as_tibble() %>%
  pull()

newdeaths_by_county <-
  deaths_df %>%
  filter(!is.na(Date)) %>%
  arrange(Date) %>%
  tail(n = 2) %>%
  select(-all_of(non_county_fields))

newdeaths_by_county <-
  (newdeaths_by_county[2, ] - newdeaths_by_county[1, ]) %>%
  as_tibble() %>%
  add_column(new_date, .before = 1) %>%
  rename(Date = new_date) %>%
  gather(Date) %>%
  rename(County = Date) %>%
  rename(newdeaths = value)

### Get the number of total deaths by county
recovered_by_county <-
  recover_df %>%
  filter(Date > "2020-03-31") %>%
  arrange(Date) %>%
  tail(n = 1) %>%
  select("Date", all_of(all_counties)) %>%
  gather(Date) %>%
  rename(County = Date) %>%
  rename(recovered = value)

### Still sick
sick_by_county <-
  totalinfected_by_county %>%
  left_join(recovered_by_county, by = "County") %>%
  left_join(totaldeaths_by_county, by = "County") %>%
  mutate(sick = totalinfected - recovered - totaldeaths) %>%
  select(County, sick)

### Combine the data
stats_df <-
  totalinfected_by_county %>%
  left_join(newinfected_by_county, by = "County") %>%
  left_join(totaldeaths_by_county, by = "County") %>%
  left_join(newdeaths_by_county, by = "County")   %>%
  left_join(sick_by_county, by = "County") %>%
  left_join(recovered_by_county, by = "County")

### Combine all the above together
counties <-
  counties %>%
  left_join(tn_map_df, by = "GEOID") %>%
  left_join(stats_df,  by = "County") %>%
  mutate(newinfected = if_else(newinfected == -1, 0, newinfected))


################################################################################
### Map of total COVID-19 cases
################################################################################
### We want to add text with the number of infected in each county to that
### county on the map.   But we *don't* want to label counties that don't yet
### have any confirmed cases.
totalinfected_label <- counties$totalinfected
totalinfected_label[totalinfected_label == 0] <- ""

### By default use black for the font, but if the value is over 1/2's of the way
### to the max, use the white font instead as it looks nicer
frac <- 0.5 * max(counties$totalinfected)
counties$textcolor <- if_else(counties$totalinfected >= frac, "white", "black")

### Map:  COVID-19 confirmed infected per county
p_map_totalinfected <-
  ggplot(counties) +
  theme_void() +
  theme(legend.title = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, aes(fill = counties$totalinfected),
          size = geom_thickness) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = totalinfected_label),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "Total Cases") +
  scale_fill_gradientn(colours = map_palette)


################################################################################
### Map of new COVID-19 cases
################################################################################
### We want to add text with the number of infected in each county to that
### county on the map.   But we *don't* want to label counties that don't yet
### have any confirmed cases.
newinfected_label <- counties$newinfected
newinfected_label[newinfected_label <= 0] <- ""

### By default use black for the font, but if the value is over 1/2's of the way
### to the max, use the white font instead as it looks nicer
frac <- 0.5 * max(counties$newinfected)
counties$textcolor <- if_else(counties$newinfected >= frac, "white", "black")

### Map:  COVID-19 confirmed new infected per county
p_map_newinfected <-
  ggplot(counties) +
  theme_void() +
  theme(legend.title = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, aes(fill = counties$newinfected),
          size = geom_thickness) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = newinfected_label),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "New Cases") +
  scale_fill_gradientn(colours = map_palette)


################################################################################
### Map of deaths by county
################################################################################
total_deaths_label <- counties$totaldeaths
total_deaths_label[total_deaths_label == 0] <- ""

### By default use black for the font, but if the value is over 1/2's of the way
### to the max, use the white font instead as it looks nicer
frac <- 0.5 * max(counties$totaldeaths)
counties$textcolor <- if_else(counties$totaldeaths >= frac, "white", "black")

p_map_totaldeaths <- ggplot(counties) +
  theme_void() +
  theme(legend.title = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, size = geom_thickness,
          aes(fill = counties$totaldeaths)) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = total_deaths_label),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "Total Deaths") +
  scale_fill_gradientn(colours = map_palette)

new_deaths_label <- counties$newdeaths
new_deaths_label[new_deaths_label == 0] <- ""

### By default use black for the font, but if the value is over 1/2's of the way
### to the max, use the white font instead as it looks nicer
frac <- 0.5 * max(counties$newdeaths)
counties$textcolor <- if_else(counties$newdeaths > frac, "white", "black")

p_map_newdeaths <- ggplot(counties) +
  theme_void() +
  theme(legend.title = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, aes(fill = counties$newdeaths),
          size = geom_thickness) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = new_deaths_label),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "New Deaths") +
  scale_fill_gradientn(colours = map_palette)


################################################################################
### Map of deaths per 100k by county
################################################################################
counties$deaths_per100k <- round(100000 * counties$totaldeaths /
                                   counties$POP2018, 3)
deaths100k_label <- round(counties$deaths_per100k, 1)
deaths100k_label[deaths100k_label == 0] <- ""

### By default use black for the font, but if the value is over 1/2's of the way
### to the max, use the white font instead as it looks nicer
frac <- 0.5 * max(counties$deaths_per100k)
counties$textcolor <- if_else(counties$deaths_per100k >= frac, "white", "black")

p_map_totaldeaths_100k <- ggplot(counties) +
  theme_void() +
  theme(legend.title = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, aes(fill = counties$deaths_per100k),
          size = geom_thickness) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = deaths100k_label),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "Total Deaths per 100k") +
  scale_fill_gradientn(colours = map_palette)


################################################################################
### Map of total COVID-19 cases per 100k
################################################################################
counties$infected_percapita <-
  round(100000 * counties$totalinfected / counties$POP2018, 0)

frac <- 0.5 * max(counties$infected_percapita)
counties$textcolor <- if_else(counties$infected_percapita >= frac,
                              "white", "black")

infected_label_per  <- round(counties$infected_percapita, 1)
infected_label_per[infected_label_per == 0] <- ""

p_map_totalinfected_per <- ggplot(counties) +
   theme_void() +
   theme(plot.title = element_text(hjust = 0.5)) +
   theme(plot.caption = element_text(hjust = 0.5)) +
   theme(legend.title = element_blank()) +
   theme(legend.position = "none") +
   geom_sf(data = counties$geometry, size = geom_thickness,
             aes(fill = counties$infected_percapita)) +
   geom_text(data = counties, size = map_fontsize, fontface = "bold",
             color = counties$textcolor,
             aes(x     = county_centers$x,
                 y     = county_centers$y,
                 label = infected_label_per),
               nudge_x = counties$nudge_x,
               nudge_y = counties$nudge_y) +
   labs(title = "Total Cases per 100k") +
     scale_fill_gradientn(colours = map_palette)


################################################################################
### Map of people still sick (not recovered) by county
################################################################################
frac <- 0.5 * max(counties$sick)
counties$textcolor <- if_else(counties$sick >= frac, "white", "black")

sick_label_per  <- counties$sick
sick_label_per[sick_label_per == 0] <- ""

p_map_sick <- ggplot(counties) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.caption = element_text(hjust = 0.5)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, size = geom_thickness,
          aes(fill = counties$sick)) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = sick_label_per),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "Still Sick") +
  scale_fill_gradientn(colours = map_palette)


################################################################################
### Map of people recovered by county
################################################################################
frac <- 0.5 * max(counties$recovered)
counties$textcolor <- if_else(counties$recovered >= frac, "white", "black")

recovered_label_per  <- counties$recovered
recovered_label_per[recovered_label_per == 0] <- ""

p_map_recovered <- ggplot(counties) +
  theme_void() +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.caption = element_text(hjust = 0.5)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_sf(data = counties$geometry, size = geom_thickness,
          aes(fill = counties$recovered)) +
  geom_text(data = counties, size = map_fontsize, fontface = "bold",
            color = counties$textcolor,
            aes(x     = county_centers$x,
                y     = county_centers$y,
                label = recovered_label_per),
            nudge_x = counties$nudge_x,
            nudge_y = counties$nudge_y) +
  labs(title = "Recovered") +
  scale_fill_gradientn(colours = map_palette)


################################################################################
### Math section for charts
################################################################################
infected_data <-
  infected_df %>%
  arrange(Date) %>%
  select("Date", all_of(graph_counties)) %>%
  mutate(Total = rowSums(select(., !starts_with("Date")))) %>%
  select(Date, Total) %>%
  rename(totalinfected = Total) %>%
  rename(date = Date) %>%
  filter(totalinfected != 0) %>%
  mutate(newinfected = c(totalinfected[1], diff(totalinfected))) %>%
  mutate(numdays = difftime(as.POSIXct(date),
                            as.POSIXct(as.Date(min(date)) - 1, tz = "UTC"),
                            units = "days")) %>%
  mutate(numdays = as.numeric(numdays))

hospital_data <-
  hospital_df %>%
  arrange(Date) %>%
  mutate(numdays = difftime(as.POSIXct(Date),
                            as.POSIXct(as.Date(min(Date)) - 1, tz = "UTC"),
                            units = "days")) %>%
  mutate(numdays = as.numeric(numdays)) %>%
  rename(total_hospitalizations = Hospitalizations) %>%
  mutate(new_hospitalizations = c(total_hospitalizations[1],
                                  diff(total_hospitalizations)))

recovered_data <-
  recover_df %>%
  rename(Total = Recovered) %>%
  select("Date", all_of(graph_counties)) %>%
  mutate(Total = rowSums(select(., !starts_with("Date")))) %>%
  select(Date, Total) %>%
  rename(totalrecovered = Total) %>%
  rename(date = Date) %>%
  filter(totalrecovered != 0) %>%
  arrange(date) %>%
  mutate(newrecovered = c(totalrecovered[1], diff(totalrecovered))) %>%
  mutate(numdays = difftime(as.POSIXct(date),
                            as.POSIXct(as.Date(min(date)) - 1, tz = "UTC"),
                            units = "days")) %>%
  mutate(numdays = as.numeric(numdays))

deaths_data <-
  deaths_df %>%
  na.omit() %>%
  arrange(Date) %>%
  select("Date", all_of(graph_counties)) %>%
  mutate(Total = rowSums(select(., !starts_with("Date")))) %>%
  mutate(New = c(1, diff(Total))) %>%
  select(Date, Total, New) %>%
  filter(Total != 0) %>%
  rename(totaldeaths = Total) %>%
  rename(newdeaths = New) %>%
  rename(date = Date) %>%
  mutate(numdays = difftime(as.POSIXct(date),
                            as.POSIXct(as.Date(min(date)) - 1, tz = "UTC"),
                            units = "days")) %>%
  mutate(numdays = as.numeric(numdays))

sick_data <-
  infected_data %>%
  left_join(recovered_data, by = "date") %>%
  left_join(deaths_data,   by = "date") %>%
  select(date, totalinfected, totaldeaths, totalrecovered) %>%
  filter(date >= as.Date("2020-04-10")) %>%
  filter(!is.na(totaldeaths)) %>%
  filter(!is.na(totalrecovered)) %>%
  filter(!is.na(totalinfected)) %>%
  mutate(current_sick = totalinfected - totaldeaths - totalrecovered) %>%
  select(date, current_sick)


################################################################################
### Graph:  COVID-19 Total Confirmed Cases
################################################################################
### Use non-linear fitting to do the exponential fit
fit_totalinfected <-
  nls(totalinfected ~ a * exp(b * numdays),
      data = infected_data,
      start = list(a = 334,
                   b = 0.0677))

### Extract the regression coefficents so we can add the regression formula
### to the graph
intercept  <- coef(fit_totalinfected)[1]
d_constant <- coef(fit_totalinfected)[2]

### Add the fit to our dataframe
infected_data <-
  infected_data %>%
  mutate(fit_totalinfected = intercept *
           exp(d_constant * infected_data$numdays))

### Compute the doubling time
doubling_time <- log(2) / d_constant

### Compute the "order-of-magnitude" time
magnitude_time <- log(10) / d_constant

### Growth rate
growth_rate <- exp(log(2) / doubling_time) - 1

### What's the current number of total infected cases
total_num <-
  infected_data %>%
  tail(n = 1) %>%
  select("totalinfected") %>%
  pull()

### What's the date of the last data point
data_date <- infected_data  %>% select("date") %>% tail(n = 1) %>% pull()

p_total <- ggplot(data = infected_data) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_point(data = infected_data, shape = 19, size = 0.5,
             aes(x = as.Date(date), y = totalinfected), color = "black") +
  geom_line(data = infected_data, color = graph_color,
            aes(x = as.Date(date), y = fit_totalinfected)) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  graph_log10_opts1 +
  graph_log10_opts2 +
  labs(title = paste("Total Cases: ", total_num, sep = ""), x = "", y = "")


################################################################################
### Graph:  COVID-19 New Confirmed Cases
################################################################################
### If we assume infected = a * B^(c * days), then
###
###  d(infected)/dt = c * log_e(B) * infected
###
### We can use the same constants from the total fit on the new fit.   If new
### cases are still growing exponentially, it will work.   And when it doesn't
### work, that's a sign that we're either "flattening the curve", or something
### else weird is going on

### Add the fit to our dataframe
infected_data <-
  infected_data %>%
  mutate(fit_newinfected = fit_totalinfected * d_constant)

### What's the current number of total infected cases
new_num <- infected_data %>% tail(n = 1) %>% select("newinfected") %>% pull()
SMA_num <- infected_data$newinfected %>% SMA(n = 7) %>% tail(n = 1) %>% round()

### Graph:  new COVID-19 Cases
p_new <- ggplot(data = infected_data) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_point(data = infected_data, shape = 19, size = 0.5,
             aes(x = as.Date(date), y = newinfected), color = "black") +
  geom_line(data = infected_data, color = graph_color,
            size = 1.0, aes(x = as.Date(date) - 3,
                            y = SMA(newinfected, n = 7))) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  graph_log10_opts1 +
  graph_log10_opts2 +
  labs(title = paste("New Cases: ", new_num, ", SMA: ", SMA_num, sep = ""),
       x = "", y = "")


################################################################################
### Graph:  COVID-19 Cases [ log(Total) vs log(New) ]
### Reference:  https://aatishb.com/covidtrends/
################################################################################
p_total_new <- ggplot(data = infected_data,
                      aes(x = totalinfected, y = newinfected)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.position = "none") +
  geom_point(data = infected_data, shape = 19, size = 0.5,
             aes(x = totalinfected, y = newinfected), color = "black") +
  geom_smooth(method = "loess", formula = y ~ x,
              se = FALSE, size = 0.5, color = graph_color) +
  scale_x_log10(breaks = c(1, 10, 100, 1000)) +
  scale_y_log10(breaks = c(1, 10, 100, 1000)) +
  annotation_logticks() +
  labs(title = "Log(New vs Total)", x = "", y = "")


################################################################################
### Graph:  COVID-19 Hospitalizations
################################################################################
hospital_num <-
  hospital_data %>%
  arrange(Date) %>%
  tail(n = 1) %>%
  select("total_hospitalizations") %>%
  pull()

p_total_hospital <- ggplot(data = hospital_data,
                           aes(x = as.Date(Date), y = total_hospitalizations)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_bar(stat = "identity", fill = graph_color) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  labs(title = paste("TN Hospitalizations: ", hospital_num, sep = ""),
       x = "", y = "")


################################################################################
### Graph:  COVID-19 Total Recovered
################################################################################\
new_num <-
  recovered_data %>%
  tail(n = 1) %>%
  select("totalrecovered") %>%
  pull()

p_recover <- ggplot(data = recovered_data,
                    aes(x = as.Date(date), y = totalrecovered)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_point(data = recovered_data, shape = 19, size = 0.5,
             aes(x = as.Date(date), y = totalrecovered), color = "black") +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  labs(title = paste("Total Recovered: ", new_num, sep = ""), x = "", y = "")


################################################################################
### Graph:  COVID-19 Total Infected by Age
################################################################################
p_age_infected <- ggplot(data = age_df, aes(x = infected_statewide, y = age)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  geom_age_orientation +
  geom_bar(stat = "identity", fill = graph_color) +
  labs(title = "TN Infected by Age", x = "", y = "")


################################################################################
### Graph:  COVID-19 Total Infected by Age
################################################################################
p_age_deaths <- ggplot(data = age_df, aes(x = deaths_statewide, y = age)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  geom_age_orientation +
  geom_bar(stat = "identity", fill = graph_color) +
  labs(title = "TN Deaths by Age", x = "", y = "")


################################################################################
### Graph:  COVID-19 Total/New Deaths
################################################################################
### Use non-linear fitting to do the exponential fit
fit_totaldeaths <- nls(totaldeaths ~ a * exp(b * numdays),
                  data = deaths_data, start = list(a = 12.51, b = 0.082))

### Extract the regression coefficents so we can add the regression formula
### to the graph
intercept  <- coef(fit_totaldeaths)[1]
d_constant <- coef(fit_totaldeaths)[2]

### Add fit for total/new deaths
deaths_data <-
  deaths_data %>%
  mutate(fit_totaldeaths = intercept * exp(d_constant * numdays)) %>%
  mutate(fit_newdeaths   = fit_totaldeaths * d_constant)

p_totaldeaths <- ggplot(data = deaths_data) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_point(shape = 19, size = 0.5,
             aes(x = as.Date(date), y = totaldeaths), color = "black") +
  geom_line(color = graph_color,
            aes(x = as.Date(date), y = fit_totaldeaths)) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  graph_log10_opts1 +
  graph_log10_opts2 +
  labs(title = paste("Total Deaths: ", max(deaths_data$totaldeaths),  sep = ""),
       x = "", y = "")

p_newdeaths <- ggplot(data = deaths_data) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_point(shape = 19, size = 0.5,
             aes(x = as.Date(date), y = newdeaths), color = "black") +
  geom_line(color = graph_color,
            aes(x = as.Date(date), y = fit_newdeaths)) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  graph_log10_opts1 +
  graph_log10_opts2 +
  labs(title = paste("New Deaths: ", max(deaths_data$newdeaths),  sep = ""),
       x = "", y = "")


################################################################################
### Graph:  COVID-19 Current Sick
################################################################################
new_num <- sick_data %>% tail(n = 1) %>% select("current_sick") %>% pull()

p_sick <- ggplot(data = sick_data, aes(x = as.Date(date), y = current_sick)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(vjust = 0.7)) +
  theme(axis.text.x = element_text(hjust = 0.8)) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  geom_point(data = sick_data, shape = 19, size = 0.5,
             aes(x = as.Date(date), y = current_sick), color = "black") +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d") +
  labs(title = paste("Currently Sick: ", new_num, sep = ""), x = "", y = "")


################################################################################
### Arrange the graphs into the final product
################################################################################
# What's the date of the last data point
data_date <- infected_df %>% arrange(Date) %>% tail(n = 1) %>% pull("Date")

# Title
title_string <- paste("COVID-19 Confirmed Cases in ", location,
                      " [", as.Date(data_date), "]", sep = "")

# Footer
footer_string <- paste("Source:  Tennessee Department of Health\n",
                       "https://www.tn.gov/health/cedep/ncov.html", sep = "")

grob_charts <-
  plot_grid(p_total,        p_new,        p_totaldeaths,    p_recover,   p_sick,
            p_age_infected, p_age_deaths, p_total_hospital, p_total_new, NA,
            ncol = 5, nrow = 2, align = "hv")

grob_maps <-
  plot_grid(p_map_recovered,    p_map_totalinfected,
            p_map_sick,         p_map_totalinfected_per,
            p_map_totaldeaths,  p_map_newinfected,
            ncol = 2, nrow = 3, align = "hv")

charts <- plot_grid(grob_maps,
                    grob_charts,
                    ncol = 1, nrow = 2,
                    rel_widths  = c(1.2, 1),
                    rel_heights = c(1, 1))

title <- ggdraw() + draw_label(title_string, fontface = "bold")
Sys.sleep(2)

footer <- ggdraw() + draw_label(footer_string, size = 10)
Sys.sleep(2)

final <- plot_grid(title, charts, footer, nrow = 3,
                   rel_heights = c(0.05, 1, 0.05))
Sys.sleep(2)

print(final)

### Save image to file and we're done!
scale <- 1.3
width <- round(1300 * scale)
height <- round(900 * scale)
png("TN_COVID_Viridis.png", width = width, height = height)
plot(final)
dev.off()