library(tidyverse)
library(sf)
library(leaflet)
library(leafpop)
library(formattable)
library(ggiraph)
# install.packages("devtools")
# devtools::install_github("UrbanInstitute/urbnmapr")
library(urbnmapr) # change to include AK/HI on map
library(cowplot) 
library(RcppRoll) 
library(htmlTable)
library(ggtext)
library(patchwork)
library(Hmisc)
library(httr)
library(rvest)
# require(RSocrata); require(plyr); require(readr); require(here)
# library(rvest)
# library(albersusa)

# Default ggplot theme
theme_set(theme_bw() +
            theme(axis.title.x = element_text(size = rel(1.25)),
                  axis.title.y = element_text(size = rel(1.25)),
                  plot.caption = element_text(face = "italic", size = 9),
                  legend.position = "bottom",
                  legend.background = element_rect(fill = "transparent"), 
                  panel.background = element_rect(fill = "transparent"),
                  panel.border = element_blank(),
                  axis.line = element_line(color="black"),
                  panel.grid.major.x = element_blank(),
                  panel.grid.minor.x = element_blank(),
                  panel.grid.minor.y = element_blank(),
                  panel.grid.major.y = element_line(color = "grey90", linetype = 2)
            ))

Last Data Update: 2023-12-19 21:31:13.102324. All source documents for this report can be found on github

This document updates weekly each Wednesday after the TN State DOH releases data (usually Wednesday at ~4PM). The CDC releases data each weekday on their site.

For official, weekly TN State Dept of Health Updates, go to the Critical Indicator Report.

VACCINE SIGN-UP

For a reliable discussion of COVID-19, vaccines, and emerging variants, listen to Dr. Aronoff’s Medicine Grand Rounds from August 2021.

All ages over 6 months are now eligible for COVID vaccination- Go to one of these sites to sign up:

Disclaimer

Update 5/10/2021 New Cases now reported only by the date of detection. Will update with details once it becomes clear how these results will be mixed with the prior reporting method (date of report).

The Tennessee state information is simply a visual summary of the data directly downloaded from TN state Dept of Health, by date reported (rather than date of test). I do not attempt to model, make predictions, or give advice on how to use this data. If you want that sort of thing, go to the VUMC Health Policy site, which has developed some very good local models and Advisory Memos for Tennessee.

For the most up-to-date numbers visit the TN state Dept of Health or local health department sites. In particular, the Metro. Nashville Dashboard has a very nice summary.

Other useful sites for following trends in Tennessee and nationally:

TN Dept of Health Data

These data are taken directly from the TN state Dept of Health. Previously the DOH reported results by date of report, rather than date of test. They transitioned to reporting results by Date of Test on 5/10/2021, although it is not clear if data reported prior to that will be cleaned up.

Hospitalizations: A current summary of Hospitalizations can be found on the TN DOH site.

# READ IN AND format data ####
tndoh_daily <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-Daily-Case-Info.xlsx"))
names(tndoh_daily) <- tolower(names(tndoh_daily))



# remove the new_test results which our bulk reported or adjusted
tndoh_daily <-
  tndoh_daily %>%  
  mutate(date = lubridate::ymd(date)) %>%  
  arrange(date) %>% 
  filter(date > lubridate::ymd("2000-01-01")) # State occasionally releases data with default SAS date of "1960-01-01"for missing values

# min(tndoh_daily$date)
# tndoh_daily <-
#   tndoh_daily %>%
#   filter(date > lubridate::ymd("2020-01-01")) # data error  introduced; filter out erroneous entry for 1960-01-01
   
tndoh_daily[tndoh_daily$date %in% lubridate::ymd(c("2020-07-13", "2020-07-11", "2020-06-12")), "new_tests"] <- NA_real_


# STATE POPULATION DATA ####
# us_data <-  readRDS(here::here("data_rds/us_data.rds"))
us_data <- get_urbn_map("states", sf = TRUE)
us_data <- 
  us_data %>% 
  mutate(name = state_name,
         stusps = state_abbv)


# state_pop <- readxl::read_excel(here::here("state_populations.xlsx"))
# state_pop$region <- gsub("\\.", "", state_pop$region)
# state_pop <-
#   state_pop %>% 
#   pivot_longer(cols = -c(region, census_2010),
#                names_to = "pop_year",
#                names_prefix = "pop_",
#                values_to = "pop_estimate")
# saveRDS(state_pop, here::here("data_rds/state_pop.rds"))

state_pop <- readRDS( here::here("data_rds/state_pop.rds"))
state_pop_2019 <-
  state_pop %>%  
  filter(pop_year == 2019) %>%  
  filter(region %in% state.name)
TN_pop2019 <- state_pop_2019$pop_estimate[state_pop_2019$region == "Tennessee"]

# vaccine state info
tndoh_vax_state <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_STATE_SUMMARY.XLSX"))
names(tndoh_vax_state) <- tolower(names(tndoh_vax_state))
tndoh_vax_state <-
  tndoh_vax_state %>% 
  mutate(date = lubridate::ymd(date),
         vaccine_count_pct = vaccine_count/TN_pop2019*100,
         recipient_count_pct = recipient_count/TN_pop2019*100,
         recip_fully_vacc_pct = recip_fully_vacc/TN_pop2019*100) 

# add trend line calculated by STL (Seasonal and Trend decomposition using Loess) 
tndoh_daily$new_cases_ts <- ts(tndoh_daily$new_cases, frequency = 7, start = 1)
tndoh_daily$new_cases_trend.stl <- stl(tndoh_daily$new_cases_ts, s.window = 7)$time.series[, "trend"]

# add Deaths trend line calculated by STL (Seasonal and Trend decomposition using Loess) 
tndoh_daily$new_deaths_ts <- ts(tndoh_daily$new_deaths_by_dod, frequency = 7, start = 1)
tndoh_daily$new_deaths_trend.stl <- stl(tndoh_daily$new_deaths_ts, s.window = 7)$time.series[, "trend"]


# create annotations
tndoh_daily <-
  left_join(tndoh_daily, tndoh_vax_state,
            by = "date") %>%  
  mutate(pos_tests_new = pos_tests -lag(pos_tests, 1), 
         pos_pct =   pos_tests/total_tests,
         pos_pct_today = pos_tests_new/new_tests,
         hospitalized_tot = replace_na(total_hosp, 0),
         hospitalized_new = total_hosp - lag(total_hosp),
         new_cases_100k = new_cases/TN_pop2019*100000,
         new_cases_7davg = RcppRoll::roll_meanr(new_cases, 7),
         new_cases_7davg_100k = new_cases_7davg/TN_pop2019*100000,
         new_deaths_7davg = RcppRoll::roll_meanr(new_deaths_by_dod, 7),
         new_hosp_7davg = RcppRoll::roll_meanr(new_hosp, 7),
         new_tests_7davg = RcppRoll::roll_meanr(new_tests, 7),
         pos_pct_7davg = RcppRoll::roll_meanr(pos_tests_new, 7)/RcppRoll::roll_meanr(new_tests, 7),
         tt_daily = sprintf("<tr><td><b>%s</b></td></tr><br>
                     <tr><td><b>New Cases: %s</b></td></tr><br>
                     <tr><td><b>New Cases/100k: %s</b></td></tr><br>
                     <tr><td><b>New Hospitalizations: %s</b></td></tr><br>
                     <tr><td><b>New Deaths: %s</b></td></tr><br>
                     <tr><td><b>Test+ %% Today: %1.1f%%</b></td></tr><br>
                     <tr><td>New Cases 7d Avg/100k: %s</td></tr><br>
                     <tr><td>New Cases trend: %s</td></tr><br>
                     <tr><td>Deaths 7d Avg: %1.1f</td></tr><br>
                     <tr><td>Test+ %% 7d Avg: %1.1f%%</td></tr><br>
                     <tr><td>Case Fatality Ratio: %1.2f%%</td></tr><br>
                     <tr><td>New Vaccinations: %s</td></tr><br>
                     <tr><td>Rec 1st Vaccination %%Pop: %1.2f%%</td></tr>",
                     paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date), ", ", lubridate::year(date), " (", lubridate::wday(date, label = T), ")"), 
                     comma(new_cases, digits = 0), 
                     comma(new_cases_100k, digits = 0), 
                     comma(new_hosp, 0), 
                     comma(new_deaths_by_dod, digits = 0), 
                     pos_pct_today*100, 
                     comma(new_cases_7davg_100k, 1),
                     comma(new_cases_trend.stl/TN_pop2019*100000, 1),
                     comma(new_deaths_7davg, digits = 0), 
                     pos_pct_7davg*100, 
                     total_deaths_by_dod/total_cases*100,
                     comma(new_vaccine_count, 0),
                     recipient_count_pct) ) 


tndoh_daily_limited <-
  tndoh_daily %>% 
  select(date, total_cases, total_hosp, total_deaths_by_dod) %>%  
  pivot_longer(cols = c(total_cases, total_hosp, total_deaths_by_dod),
               names_to = "measure",
               values_to = "count") %>%  
  mutate(var.f = factor(measure, levels = c("total_cases", "total_hosp", "total_deaths_by_dod"),
                        labels = c("Cases", "Hospitalizations", "Deaths")))  
title_tnsdc <- cowplot::ggdraw() + 
  cowplot::draw_label(paste0("COVID-19 in Tennessee (as of ", max(tndoh_daily$date), ")") , fontface='bold')

tndoh_daily_long <-
  tndoh_daily %>%  
  select(date, new_cases, new_deaths_by_dod, new_hosp) %>%   
  pivot_longer(cols = c(new_cases, new_deaths_by_dod, new_hosp),
               names_to = "var", 
               values_to = "value") %>%  
  mutate(var.f = factor(var, levels = c("new_cases", "new_deaths_by_dod", "new_hosp"), 
                        labels = c("Cases", "Hospitalizations", "Deaths")))

label_pos_pct <-
    with(tndoh_daily,
         sprintf("Overall Positive: %1.1f%%\nToday's Positive: %1.1f%%", 
                 pos_pct[which(date == max(date))]*100,
                 pos_pct_today[which(date == max(date))]*100))



# AGE DATA ####
tndoh_age <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-Age.xlsx"))
names(tndoh_age) <-  tolower(names(tndoh_age))
tndoh_age <-
  tndoh_age %>% 
  filter(age_range != "Pending") %>%  
  mutate(date = lubridate::ymd(date),
         age_range = gsub(" years", "", age_range),
         age.f = factor(age_range),
         deaths_case_rate = case_when(age.f == "Pending" ~ NA_real_,
                                      TRUE ~ ar_totaldeaths/ar_totalcases*100))

# TN age groups 2019 estimates(taken from https://www.census.gov/data/tables/time-series/demo/popest/2010s-state-detail.html)
tn_age_groups <-
  tibble::tibble(age.f = factor(c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"),
                                    levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
                 populationtn_byage_2019 = c(908244, 858307, 943700, 879850, 854356, 893802, 782810, 481538, 226567))
                                         # c(821984, 855582, 942211, 887898, 850125, 893578, 808164, 516865, 252767)) # OLD ESTIMATES; RECALCULATED
tn_age_groups2 <-
  tibble::tibble(age.f2 = factor(c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"), # 16-20 needs to match other age category
                        levels = c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
                 populationtn_byage_2019_B = c(408605, 585900, 343637, 428409, 943700, 879850, 854356, 893802, 782810, 481538, 226567))


tndoh_age <- 
  left_join(tndoh_age,
          tn_age_groups,
          by = "age.f") %>%  
  group_by(age.f) %>%  
  mutate(ar_new_cases_7davg = round(RcppRoll::roll_meanr(ar_newcases, 7), 1),
         ar_new_cases_100k = ar_newcases/populationtn_byage_2019*100000,
         ar_new_cases_100k_7davg = round(RcppRoll::roll_meanr(ar_newcases, 7)/populationtn_byage_2019*100000, 1),
         ar_totaldeaths_100k    = ar_totaldeaths/populationtn_byage_2019*100000,
         ar_new_deaths_7davg = round(RcppRoll::roll_meanr(ar_newdeaths, 7), 1),
         ar_new_deaths_100k = ar_newdeaths/populationtn_byage_2019*100000,
         ar_new_deaths_100k_7davg = round(RcppRoll::roll_meanr(ar_newdeaths, 7)/populationtn_byage_2019*100000, 1)) %>%
  ungroup()

tndoh_age <- 
  tndoh_age %>%
  group_by(date) %>% 
  mutate(total_new_cases_100k_7davg = sum(ar_new_cases_100k_7davg, na.rm = T),
         ar_proportion_100k_7davg = ar_new_cases_100k_7davg/total_new_cases_100k_7davg,
         ar_newcases_sum_7davg = sum(ar_new_cases_7davg, na.rm = T),
         ar_proportion_7davg = ar_new_cases_7davg/ar_newcases_sum_7davg,
         ar_deaths_sum_7davg = sum(ar_new_deaths_7davg, na.rm = T),
         ar_newdeaths_proportion_7davg = ar_new_deaths_7davg/ar_deaths_sum_7davg) %>% 
  ungroup()


# VACCINE by age
tndoh_vax_age <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_AGE_GROUPS.XLSX"),
                                      col_types = c("text", "date", rep("numeric", 8)))
names(tndoh_vax_age) <- tolower(names(tndoh_vax_age))
tndoh_vax_age <-
  tndoh_vax_age %>%
  filter(!grepl("pending", age_group, ignore.case = T )) %>%
  mutate(date = lubridate::ymd(date),
         age.f2 = factor(age_group,
                        levels = c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"),
                        labels = c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"))) 
# tndoh_vax_age <-
#   tndoh_vax_age %>%
#   filter(!grepl("pending", age_group, ignore.case = T )) %>%  
#   mutate(date = lubridate::ymd(date),
#          age.f = factor(age_group,
#                         levels = c("0-10", "11-20", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"), # 16-20 needs to match other age category
#                         labels = c("0-10", "11-20", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")))
tndoh_vax_age <-
  left_join(tndoh_vax_age, tn_age_groups2,
            by = "age.f2") %>%
  rowwise() %>%
  mutate(ar_totalvax_pct = vaccine_count/populationtn_byage_2019_B*100,
         ar_recipientvax_pct = recipient_count/populationtn_byage_2019_B*100) %>%
  ungroup() 

tndoh_age_last <- # FOR AGE PROPORTION PLOT
  tndoh_age %>% 
  filter(!is.na(ar_proportion_7davg)) %>% 
  filter(date == max(date)) %>% 
  arrange(desc(age_range)) %>% 
  mutate(sec_axis.ticks_cases = cumsum(ar_proportion_7davg),
         sec_axis.label_cases = sprintf("<b>%s:</b> %1.1f%%", age_range, ar_proportion_7davg*100),
         sec_axis.ticks_deaths = cumsum(ar_newdeaths_proportion_7davg),
         sec_axis.label_deaths = sprintf("<b>%s:</b> %1.1f%%", age_range, ar_newdeaths_proportion_7davg*100)) %>% 
  select(age_range, ar_proportion_7davg, ar_newdeaths_proportion_7davg,
         sec_axis.ticks_cases, sec_axis.label_cases, 
         sec_axis.ticks_deaths, sec_axis.label_deaths) %>% 
  mutate(sec_axis.ticks_cases_mid = (sec_axis.ticks_cases + lag(sec_axis.ticks_cases, default = 0))/2,
         sec_axis.ticks_deaths_mid = (sec_axis.ticks_deaths + lag(sec_axis.ticks_deaths, default = 0))/2)


# most recent vaccination by age
vax_by_age_mostrecent <-
  tndoh_vax_age %>% 
  filter(!is.na(vaccine_count)) %>% 
  filter(date == max(date, na.rm = T)) %>% 
    select(date,age.f2, recipient_count, populationtn_byage_2019_B, ar_recipientvax_pct) 
mostrecent_vax_agedate <- max(vax_by_age_mostrecent$date)  

vax_by_age_mostrecent_refactor <- select(vax_by_age_mostrecent, -populationtn_byage_2019_B)
vax_by_age_mostrecent_refactor <- vax_by_age_mostrecent_refactor[!vax_by_age_mostrecent_refactor$age.f2 %in% c("5-11", "12-15", "16-20"), ]
vax_by_age_mostrecent_refactor <-
  vax_by_age_mostrecent_refactor %>% 
  mutate(age.f = factor(age.f2,
                        levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"))) %>% 
  bind_rows(
    tibble(age.f = factor(c("0-10", "11-20"),
                          levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
           recipient_count = c(vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "5-11"]*5/6, # 0 for <5 not here
                               vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "5-11"]*1/6 +
                                 vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "12-15"] +
                                 vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "16-20"]))) %>% 
  left_join(tn_age_groups, by = "age.f")    %>% 
  mutate(ar_recipientvax_pct = recipient_count/populationtn_byage_2019*100)

# Deaths by age table ####
age_death_table <-
  tndoh_age %>%
  filter(date == max(date, na.rm = T)) %>% 
  select(date, age.f,populationtn_byage_2019, ar_totalcases, ar_totaldeaths, ar_totaldeaths_100k) %>% 
  mutate(#ar_totalcases = comma(ar_totalcases, 0),
    ar_totaldeaths_pct = sprintf("%1.2f%%", ar_totaldeaths_100k/1000),
    ar_totaldeaths_100k = sprintf("%1.2f", ar_totaldeaths_100k)) %>% 
  left_join({select(vax_by_age_mostrecent_refactor, 
       age.f, recipient_count, ar_recipientvax_pct )},
       by = "age.f")

# COUNTY DATA ####
# tndoh_countynew <- readxl::read_excel(here::here("tn_data_center/TN_COVID19_CountyDaily.xlsx"))
TN_counties <- readRDS(here::here("data_rds/TN_counties.rds")) # see below for source

tndoh_county_coln <-   length(names(readxl::read_excel(here::here("data_tndoh/Public-Dataset-County-New.xlsx"), n_max = 2)))

tndoh_countynew <-  readxl::read_excel(here::here("data_tndoh/Public-Dataset-County-New.xlsx"),
                                       col_types = c("text","date", rep("numeric", tndoh_county_coln-2))) 

names(tndoh_countynew) <-  tolower(names(tndoh_countynew))

# create STL new case trend
tndoh_countynew <-
  tndoh_countynew %>% 
  group_by(county) %>% 
  arrange(date) %>% 
  nest() %>% 
  mutate(new_cases_county_ts = map(data, ~ts(.x$new_cases, frequency = 7, start = 1)),
         new_cases_county_trend = map(new_cases_county_ts, ~stl(.x, s.window = 7)$time.series[, "trend"])) %>% 
  select(county, data, new_cases_county_trend) %>% 
  unnest(cols = c(data, new_cases_county_trend)) 


tndoh_countynew <-
  tndoh_countynew %>%  
  mutate(date = lubridate::ymd(date),
         pct_pos = total_cases/(total_tests)*100,
         county_lower = tolower(county)) %>%  
  group_by(county) %>% 
  arrange(date) %>%  
  mutate(new_cases_county_4davg = round(RcppRoll::roll_meanr(new_cases, 4), 1),
         new_cases_county_7davg = round(RcppRoll::roll_meanr(new_cases, 7), 1),
         max_new_cases_trend = max(new_cases_county_trend, na.rm = T),
         new_cases_county_trend_ratio = new_cases_county_trend/max_new_cases_trend,
         new_deaths_county_7davg = round(RcppRoll::roll_meanr(new_deaths, 7), 1),
         new_testpos_pct   = round(new_pos_tests/new_tests*100, 1),
         new_testpos_pct_7davg   = round(RcppRoll::roll_meanr(new_pos_tests, 7)/RcppRoll::roll_meanr(new_tests, 7)*100, 1),
         new_hospitalized_county_7davg = round(RcppRoll::roll_meanr(new_hospitalized, 7), 1),
         new_cases_county_7davg_1wkdelta = (new_cases_county_7davg - lag(new_cases_county_7davg, 7))/lag(new_cases_county_7davg, 7) 
  ) %>%  
  ungroup()

# TN_data <- tigris::counties(state = "TN", cb = TRUE)
# names(TN_data) <- tolower(names(TN_data))
# TN_counties <- sf::st_as_sf(TN_data) %>%  
#     mutate(county_lower = tolower(as.character(name)))

# TN COUNTY POPULATION FROM US CENSUS 2019 ESTIMATES
# tn_county_pop <- read_csv(here::here("data_csv/TN 2019 county pop estimates.csv"))
# tn_county_pop <-
#   tn_county_pop %>%  
#   mutate(county_lower = tolower(county)) %>%  
#   select(-county)
# TN_counties <-
#   left_join(TN_counties,
#             tn_county_pop,
#             by = "county_lower")


# saveRDS(TN_counties, here::here("data_rds/TN_counties.rds"))
# TN_counties <- readRDS(here::here("data_rds/TN_counties.rds"))

tndoh_countynew <-
  left_join(tndoh_countynew,
            TN_counties,
            by = "county_lower") %>%  
  mutate(total_cases_per100kpop = total_cases/population*100000,
         new_cases_per100kpop_7davg   = new_cases_county_7davg /population*100000,
         new_cases_county_trend_per100kpop   = new_cases_county_trend/population*100000,
         total_tests_pctpop = total_tests/population*100,
         new_tests_per100kpop   = new_tests/population*100000,
         total_deaths_per100kpop = total_deaths/population*100000) 

tndoh_countynew_calc <-
  tndoh_countynew %>%  
  group_by(date) %>%  
  filter(tolower(county) %in% TN_counties$county_lower) %>%  # exclude out of state and pending
  mutate(#new_testpos_pct = new_pos_tests/new_tests*100,
    new_cases_county_7davg_rank = rank(-new_cases_county_7davg),
    new_cases_per100kpop_7davg_rank = rank(-new_cases_per100kpop_7davg),
    new_cases_rank = rank(-new_cases)) %>%  
  ungroup() %>%  
  select(county, date, new_cases_county_7davg_rank, 
         new_cases_per100kpop_7davg_rank, new_cases_rank)

tndoh_countynew <-   
  left_join(tndoh_countynew,
            tndoh_countynew_calc,
            by = c("date", "county"))  


# DEMOGRAPHICS
tndoh_demographics <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-RaceEthSex.xlsx"))
names(tndoh_demographics) <-  tolower(names(tndoh_demographics))
tndoh_demographics$cat_detail <- gsub("/ ", "/", tndoh_demographics$cat_detail)
tndoh_demographics <-
  tndoh_demographics %>%  
  filter(category == "RACE",
         date == max(date)) %>%  
  dplyr::rename(race = cat_detail,
         cases_byrace = cat_totalcases) %>%  
  mutate(date = lubridate::ymd(date),
         race = factor(race,
                       levels = c("White", "Black or African American", "Asian",  "Other/Multiracial", "Pending",  "American Indian or Alaska Native", "Native Hawaiian or Other Pacific Islander"),
                       labels = c("White", "Black", "Asian",  "Other", "Pending",  "Am.Indian/Alaskan", "Hawaiian/Pacific")))

# US AND STATE DATA ####
# us_data <- tigris::states(cb = TRUE)
# us_data <-   sf::st_as_sf(us_data)
# names(us_data) <- tolower(names(us_data))
# us_data <-
#   us_data %>%
#   filter(stusps %in% state.abb)
# saveRDS(us_data, here::here("data_rds/us_data.rds"))
# us_data <-  readRDS(here::here("data_rds/us_data.rds"))
# 
# 
# 
# # state_pop <- readxl::read_excel(here::here("state_populations.xlsx"))
# # state_pop$region <- gsub("\\.", "", state_pop$region)
# # state_pop <-
# #   state_pop %>% 
# #   pivot_longer(cols = -c(region, census_2010),
# #                names_to = "pop_year",
# #                names_prefix = "pop_",
# #                values_to = "pop_estimate")
# # saveRDS(state_pop, here::here("data_rds/state_pop.rds"))
# 
# state_pop <- readRDS( here::here("data_rds/state_pop.rds"))
# state_pop_2019 <-
#   state_pop %>%  
#   filter(pop_year == 2019) %>%  
#   filter(region %in% state.name)


# US DATA FROM THE COVIDTRACKING PROJECT ####
# COVIDTRACKING PROJECT QUIT TRACKING DATA ON MARCH 7, 2021; 
# DATA SOURCE CHANGED TO CDD; TEST POSITIVITY RATES NOT AVAILABLE THERE
# covus_currentdata <- readr::read_csv(here::here("data_covidtracking/daily.csv"))
# covus_currentdata <-
#   covus_currentdata %>%  
#   filter(state %in% state.abb,
#          !state %in% c("AK", "HI") ) %>%  
#   mutate(date = lubridate::ymd(date))

# covus_currentdata <- readRDS(here::here("data_covidtracking/covus_current.rds"))
# covus_currentdata <-
#   covus_currentdata %>%  
#   left_join(us_data,
#             by = c("state" = "stusps")) %>%  
#     left_join(state_pop_2019, by = c("name" = "region")) %>%  
#   mutate(test_density_pct = total/pop_estimate*100,
#          case_density_pct = positive/pop_estimate*100)  
# 
# covus_currentdata <-
#   covus_currentdata %>% 
#   group_by(state) %>%  
#   arrange(date) %>%  
#   mutate(new_cases_state_7davg = round(RcppRoll::roll_meanr(positiveIncrease, 7), 1),
#          new_deaths_state_7davg = round(RcppRoll::roll_meanr(deathIncrease, 7), 1),
#          new_deaths_state_7davg_100k = new_deaths_state_7davg/pop_estimate*100000,
#          new_testpos_pct   = round(positiveIncrease/totalTestResultsIncrease*100, 1),
#          new_testpos_pct_7davg   = round(RcppRoll::roll_meanr(positiveIncrease, 7)/RcppRoll::roll_meanr(totalTestResultsIncrease, 7)*100, 1),
#          new_cases_state_7davg_100k = new_cases_state_7davg/pop_estimate*100000
#   ) %>%  
#   ungroup()
# 
# covus_currentdata <-
#   covus_currentdata %>% 
#   group_by(date) %>%  
#   mutate(rank_newcases7davg = rank(-new_cases_state_7davg_100k)) %>%  
#   ungroup()
# 
# state_tbl <- 
#   tibble(state_name = state.name, 
#          state = state.abb)
# 
# covus_currentdata_2163 <-
#   covus_currentdata %>%  
#   mutate(geometry = st_transform(geometry, 2163)) %>%  
#     left_join(state_tbl, by = "state") %>%  
#   mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr>
#                      <tr><td>New Cases/100k (7d Avg): %s</td></tr>
#                      <tr><td>Test %% Positive (7d Avg): %1.1f%%</td></tr>
#                      <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr>
#                      <tr><td>Total Cases: %s</td></tr>
#                      <tr><td>Total Deaths: %s</td></tr>
#                      <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr>
#                      <tr><td>CFR: %1.2f%%</td></tr>
#                      <tr><td>Case Density: %1.1f%% of pop.</td></tr>
#                      <tr><td>Testing Density: %1.1f%% of pop.</td></tr>", 
#                      toupper(state_name),  paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
#                      comma(new_cases_state_7davg_100k, digits = 1), new_testpos_pct_7davg,  rank_newcases7davg, comma(positive, digits = 0),
#                      comma(death, digits = 0), new_deaths_state_7davg_100k, death/positive*100, case_density_pct, test_density_pct)) 

# US DATA FROM THE CDC ####
# NOTE: COVIDTRACKING PROJECT QUIT TRACKING DATA ON MARCH 7, 2021; 
# TEST POSITIVITY RATES NOT AVAILABLE THERE
# CDC DATA

# cdc data download not available or moved and I can't find it 7/2023
# cdc_currentdata <- readr::read_csv(here::here("data_cdc/cdc_state_info.csv"))
# names(state.name) <- state.abb
# cdc_currentdata <-
#   cdc_currentdata %>% 
#   filter(state %in% state.abb) %>% 
#          # !state %in% c("AK", "HI") ) %>% 
#   mutate(date = lubridate::mdy(submission_date),
#          region = state.name[state])
# 
# 
# cdc_currentdata <-
#   cdc_currentdata %>%  
#   # left_join(us_data,
#   #           by = c("state" = "stusps")) %>%  
#   left_join(state_pop_2019, by = "region") 
# 
# 
# cdc_currentdata <-
#   cdc_currentdata %>% 
#   group_by(state) %>%  
#   arrange(date) %>%  
#   mutate(case_density_pct = tot_cases/pop_estimate*100,
#          new_cases_state_7davg = round(RcppRoll::roll_meanr(new_case, 7), 1),
#          new_deaths_state_7davg = round(RcppRoll::roll_meanr(new_death, 7), 1),
#          new_deaths_state_7davg_100k = new_deaths_state_7davg/pop_estimate*100000,
#          new_cases_state_7davg_100k = new_cases_state_7davg/pop_estimate*100000
#   ) %>%  
#   ungroup()
# 
# cdc_currentdata <-
#   cdc_currentdata %>% 
#   group_by(date) %>%  
#   mutate(rank_newcases7davg = rank(-new_cases_state_7davg_100k)) %>%  
#   ungroup()
# 
# 
# 
# state_tbl <-
#   tibble(state_name = state.name,
#          state = state.abb)
# 
# 
# cdc_currentdata_2163 <-
#   cdc_currentdata %>%  
#   filter(date == max(date)) %>% 
#   left_join(us_data, by = c("state" = "stusps"))  %>% 
#   mutate(geometry = st_transform(geometry, 2163)) %>%  
#   left_join(state_tbl, by = c("state", "state_name")) %>% 
#   mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr><br>
#                      <tr><td>New Cases/100k (7d Avg): %s</td></tr><br>
#                      <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr><br>
#                      <tr><td>Total Cases: %s</td></tr><br>
#                      <tr><td>Total Deaths: %s</td></tr><br>
#                      <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr><br>
#                      <tr><td>CFR: %1.2f%%</td></tr><br>
#                      <tr><td>Case Density: %1.1f%% of pop.</td></tr>",
#                      toupper(state_name),  paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
#                      comma(new_cases_state_7davg_100k, digits = 1), rank_newcases7davg, comma(tot_cases, digits = 0),
#                      comma(tot_death, digits = 0), new_deaths_state_7davg_100k, tot_death/tot_cases*100, case_density_pct)) 


  # mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr>
  #                    <tr><td>New Cases/100k (7d Avg): %s</td></tr>
  #                    <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr>
  #                    <tr><td>Total Cases: %s</td></tr>
  #                    <tr><td>Total Deaths: %s</td></tr>
  #                    <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr>
  #                    <tr><td>CFR: %1.2f%%</td></tr>
  #                    <tr><td>Case Density: %1.1f%% of pop.</td></tr>",
  #                    toupper(state_name),  paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
  #                    comma(new_cases_state_7davg_100k, digits = 1), rank_newcases7davg, comma(tot_cases, digits = 0),
  #                    comma(tot_death, digits = 0), new_deaths_state_7davg_100k, tot_death/tot_cases*100, case_density_pct)) 

# county
tndoh_vax_county <-
  readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_COUNTY_SUMMARY.XLSX"),
                       col_types = c("date", "text", rep("numeric", 8)))
    

names(tndoh_vax_county) <- tolower(names(tndoh_vax_county))
tndoh_vax_county <-
  tndoh_vax_county %>% 
  mutate(date = lubridate::ymd(date),
         county_lower = tolower(county)) %>% 
    left_join(TN_counties,
              by = "county_lower") %>% 
    rowwise() %>%  
    mutate(vaccine_count_pct = vaccine_count/population*100,
           recipient_count_pct = recipient_count/population*100) %>% 
    ungroup() %>% 
    filter(!is.na(statefp)) 
latest_vax_date <- max(tndoh_vax_county$date, na.rm = T)
tndoh_vax_county_current <- 
  left_join(
  {tndoh_vax_county %>% 
      filter(date == latest_vax_date,
             !is.na(statefp))},
  {tndoh_countynew %>%    
      filter(date == max(date)) %>% 
      select(county_lower, total_cases_per100kpop, new_cases_per100kpop_7davg)}
) %>% 
  mutate(immune_pct = recipient_count_pct + total_cases_per100kpop/1000)  

# CDC WASTEWATER DATA (added Oct 4, 2022) ----
nwss_dat <- RSocrata::read.socrata(url = "https://data.cdc.gov/resource/g653-rqe2.csv")
# nwss_dat <- content(httr::GET("https://data.cdc.gov/resource/g653-rqe2.csv")) # get data
wwtp_id_details_tn <- readRDS(here::here("data_rds/wwtp_id_details_tn.rds"))

nwss_dat <-
  nwss_dat %>% 
  mutate(original_kpid = key_plot_id) %>% 
  extract(key_plot_id, 
          into = c("nwss_source", "state", "wwtp_id", "sample_location", "sample_location_specify", "sample_type"),
          regex = "(CDC_BIOBOT|CDC_LUMINULTRA|NWSS)_([a-z]+)_([0-9]+)_(Before treatment plant|Treatment plant)_([0-9]+)*_*(.*)") %>% 
  rename(key_plot_id = original_kpid)  %>% 
  mutate(wwtp_id = as.numeric(wwtp_id), # required for merging
         date = lubridate::ymd(date))  

nwss_dat_tn <-
  nwss_dat %>% 
  filter(state == "tn") %>%
  left_join(wwtp_id_details_tn,  by = "wwtp_id")

# CDC VARIANT DATA
## ├ Get Data ----
# data source: https://covid.cdc.gov/covid-data-tracker/#variant-proportions
# download.file("https://data.cdc.gov/api/views/jr58-6ysp/rows.csv?accessType=DOWNLOAD",
#               destfile = here::here("data_cdc/cdc_variant_nowcast.csv")) # contains NOW-CAST results- all US regions
cdc_variants_data <- readr::read_csv(here::here("data_cdc/cdc_variant_nowcast.csv"))

cdc_variants_data$usa_or_hhsregion <- factor(cdc_variants_data$usa_or_hhsregion,
                                        levels = c("USA", 1:10),
                                        labels = c("USA", paste("Region", 1:10)))
cdc_variants_data$variant <- factor(cdc_variants_data$variant)
cdc_variants_data$week_ending <- lubridate::date(lubridate::mdy_hms(cdc_variants_data$week_ending))
cdc_variants_data$published_date <- lubridate::date(lubridate::mdy_hms(cdc_variants_data$published_date))

Recent Trend Update

Interactive Figure: mouse/hover over for details.

pl_new_tnsdc <-
  tndoh_daily %>% 
    filter(date >=  max(date) - 120) %>% 
  pivot_longer(cols = c(new_cases, new_hosp, new_deaths_by_dod),
               names_to = "var", 
               values_to = "value") %>%  
  # filter(date > lubridate::ymd("2021-01-01")) %>% 
  ggplot(aes(x = date)) +
  geom_rect(aes(xmin = max(tndoh_daily$date)-10, xmax = max(tndoh_daily$date) + 1, ymin = -Inf, ymax = Inf),
            fill = "grey", alpha = 0.2) +
  ggiraph::geom_col_interactive(aes(y = value, fill = var, tooltip = tt_daily, data_id = date), alpha = 0.5,
                               position = "identity") +
  geom_line(aes(y = new_cases_trend.stl),color = "red", size = 1) +
  geom_line(aes(y = new_cases_7davg), color = "red",  size = 0.5, linetype=2) +
  labs(x=NULL,
       y="Results by day reported",
       subtitle = "Results by Day, New (Total)",
       caption = "Solid line is time series trend. Dotted line is 7d avg\nResults within the past 7-10 days (shaded period) may not be fully reported",
       fill = NULL) +
  scale_fill_manual(values = c("steelblue",
                               "red",
                               "orange"),
                    labels = c(paste0("Cases: ", comma(tndoh_daily$new_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                      " (", comma(tndoh_daily$total_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
                               paste0("Deaths: ", comma(tndoh_daily$new_deaths_by_dod[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                      " (", comma(tndoh_daily$total_deaths[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
                               paste0("Hospitalizations: ", comma(tndoh_daily$new_hosp[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                      " (", comma(tndoh_daily$total_hosp[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")")
                    )) +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = sec_axis( trans = ~./TN_pop2019*10^5, name="Cases/100k Population")) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b" ) +
  theme(
    legend.position = c(.01, 1), 
    legend.justification = c(0, 1),
    legend.text = element_text(size = 10))

hover_css <- "background-color:orange;stroke:black;fill:orange;"
tooltip_css <- " background-color:rgba(255, 255, 255, 0.9);color:black;stroke:black;opacity:0.1;border:2px solid black;border-radius: 8px;font-family:Arial;font-size: 0.875em;"

ggiraph::girafe(ggobj = pl_new_tnsdc,
                width_svg = 8,
                options = list(
                  offx = 40, offy = 40,
                  opts_hover(css = hover_css),
                  opts_tooltip(css = tooltip_css)
                ) )

Deaths plotted by date-of-death. Reporting lags, so deaths within last 1-3 months may not be listed.

Limited to Last 4 months for better detail. New case trend statistic and line estimated by Seasonal and Trend decomposition using Loess (STL) decomposition to account for daily variation- results are similar to a 7-day average but better account for the cyclical variation.

Today’s Overview

# PLOTS ####
pl_tests_new <-
  tndoh_daily %>% 
  # filter(date - 90) %>% 
  ggplot(aes(x=date)) +
  geom_col(aes(y = new_tests, fill = "Tests"), alpha = 0.5) +
  geom_col(aes(y = new_cases, fill = "Cases"), alpha = 0.7) +
  labs(x=NULL,
       fill = NULL,
       y="Tests by day reported",
       subtitle = "Tests by Day, New (Total)") +
  scale_fill_manual(values = c("Cases" = "steelblue",
                               "Tests" = "grey"),
                      labels = c(paste0("Cases:", comma(tndoh_daily$new_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
                                        " (", comma(tndoh_daily$total_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
                                 paste0("Tests:", comma(tndoh_daily$new_tests[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), 
                                        " (", comma(tndoh_daily$total_tests[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"))) +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="per 100k Population")) +
  scale_x_date(date_breaks = "3 months", date_labels = "%b-%Y" ) +
  annotate(geom = "text", x=min(tndoh_daily$date)+1, y = 0.65*max(tndoh_daily$new_tests, na.rm = T), hjust = 0, 
           label = label_pos_pct, size = 3.5) +
  theme(#axis.text.x =  element_text(angle = 45, hjust = 1),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = c(.01, 1), 
        legend.justification = c(0, 1),
        legend.text = element_text(size = 10))


# age ####

pl_age <-
  tndoh_age %>% 
  filter(date == max(date)) %>% 
  ggplot() +
  geom_blank(data = tndoh_age[tndoh_age$ar_totalcases == max(tndoh_age$ar_totalcases), ], 
             aes(x=age.f, y=ar_totalcases*1.1)) +
  geom_bar(aes(age.f, ar_totalcases), alpha = 0.5, fill = "steelblue", stat = "identity", ) +
  geom_bar(aes(age.f, ar_totaldeaths), alpha = 0.5, fill = "red", stat = "identity", ) +
  geom_text(aes(age.f, ar_totalcases*1.02, vjust =0, label = comma(ar_totalcases, digits = 0)),
            color = "steelblue", size = 3) +
    geom_text(aes(age.f, ar_totaldeaths*1.01, vjust =0,  label = comma(ar_totaldeaths, digits = 0)),
              color = "red", size = 3) +
    labs(x=NULL,
       y="Total Case Count",
       subtitle = "Total Cases by Age") +
  scale_y_continuous(expand = c(0,0)) +
  theme(axis.text.x =  element_text(angle = 45, hjust = 1)) 

# multi-axis plot
pl_age_death <-
  tndoh_age %>%  
    filter(date == max(date)) %>%  
  ggplot(aes(age.f, deaths_case_rate, group = date)) +
# pl_age_death <-
#   pl_age_death +
  geom_point(color = "red") +
  geom_line(color = "red")

# Dual axis 
ggp <- ggplot_build(pl_age)
ggp2 <- ggplot_build(pl_age_death)

ylim.prim <- ggp$layout$panel_params[[1]]$y.range   # in this example, precipitation
ylim.sec <- ggp2$layout$panel_params[[1]]$y.range    # in this example, temperatureb <- diff(ylim.prim)/diff(ylim.sec)
b <- diff(ylim.prim)/diff(ylim.sec)
a <- b*(ylim.prim[1] - ylim.sec[1])

pl_age_death_final <-
  tndoh_age %>% 
  filter(date == max(date)) %>%  
  ggplot(aes(x=age.f)) +
  geom_blank(data = tndoh_age[tndoh_age$ar_totalcases == max(tndoh_age$ar_totalcases), ], 
             aes(x=age.f, y=ar_totalcases*1.1)) +
  geom_bar(aes(age.f, ar_totalcases), alpha = 0.5, fill = "steelblue", stat = "identity", ) +
  geom_bar(aes(age.f, ar_totaldeaths), fill = "red", stat = "identity", ) +
  geom_line(aes(y = a + deaths_case_rate*b, group = date), color = "red", size = 1.5) +
  geom_point(aes(y = a + deaths_case_rate*b, group = date), color = "red", size = 2, shape = 21, fill = "white") +
  geom_text(aes(age.f, ar_totalcases*1.02, vjust =0, label = comma(ar_totalcases, digits = 0)),
            color = "steelblue", size = 3) +
    geom_text(aes(age.f, ar_totaldeaths*1.01, vjust =0,  label = comma(ar_totaldeaths, digits = 0)),
              color = "red", size = 3) +
    labs(x=NULL,
       y="Total Case Count",
       subtitle = "Total Cases by Age") +
  scale_y_continuous(expand = c(0,0),
                     sec.axis = sec_axis(~ (. - a)/b, name = "Crude Case Fatality (%)")) +
  theme(axis.text.x =  element_text(angle = 45, hjust = 1),
        axis.title.y.right = element_text(color = "red"),
        axis.text.y.right = element_text(color = "red")) 


race_total <- sum(tndoh_demographics$cases_byrace[tndoh_demographics$race != "Pending"])

# AGE HEATMAP
max_newcaserate_byage <- max(tndoh_age$ar_new_cases_100k_7davg, na.rm = T)
ar_newcase_heatmap_100k <-
  tndoh_age %>%  
  filter(date >= lubridate::ymd("2020-04-01")) %>%  
  ggplot(aes(x = date,
             y = age.f)) +
  geom_tile(aes(fill = ar_new_cases_100k_7davg,
                color = ar_new_cases_100k_7davg)) +
  viridis::scale_fill_viridis(labels = scales::comma,
                              limits = c(0, max_newcaserate_byage), oob = scales::squish,
                              option = "plasma", na.value="black") +
  viridis::scale_color_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate_byage), oob = scales::squish,
                               option = "plasma", na.value="black") +
  labs(fill = "New Cases/100k/d\n(7dAvg)",
       color  = "New Cases/100k/d\n(7dAvg)",
       x = NULL, y = NULL,
       subtitle = "New Case Rate, by Age (population adjusted)") +
  scale_x_date(date_breaks = "3 months", date_labels = "%b-%Y", expand = c(0, 0)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 12),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.background = element_blank(), 
        panel.border = element_blank(),
        legend.title = element_text(size = 10),
        legend.text =  element_text(size = 8))

# fig 1 ####  
pl_TOP <- cowplot::plot_grid(pl_new_tnsdc, pl_tests_new, 
                             pl_age_death_final, ar_newcase_heatmap_100k,
                             align = "h", axis = "bt",
                             nrow = 2)
cowplot::plot_grid(title_tnsdc, pl_TOP, nrow = 2, rel_heights = c(.1, 1))

# ggsave(here::here("figures/pl_tnmap.jpg"), dpi = 600, width = 12, height = 4, units = "in")

SUMMARY BY AGE In Tennessee, 4.22% of residents >80 and 1.76% of those age 71-80 have died of COVID-19, to date.

Age Group Population (n) Cases Cases
(% of Total)
Deaths CFR% Deaths
(% of Total)
Deaths/100k Deaths (% in Age) Given 1st Vaccine (% in Age)
0-10 908244 241111 9% 18 0.007 0% 1.98 0.00% 10.9%
11-20 858307 346009 13% 28 0.008 0% 3.26 0.00% 47.3%
21-30 943700 440444 17% 222 0.050 1% 23.52 0.02% 58.2%
31-40 879850 401487 15% 599 0.149 2% 68.08 0.07% 65.1%
41-50 854356 364574 14% 1528 0.419 5% 178.85 0.18% 68.9%
51-60 893802 346939 13% 3558 1.026 12% 398.07 0.40% 75.8%
61-70 782810 260806 10% 6207 2.380 21% 792.91 0.79% 87.2%
71-80 481538 163162 6% 8465 5.188 28% 1757.91 1.76% 93.9%
81+ 226567 80898 3% 9564 11.822 32% 4221.27 4.22% 84.7%
Population within each age group taken from US Census estimates for TN, 2019
CFR=Case Fatality Rate (%)
Vaccination data last updated on 2022-08-22.

NOTES:

  • “new_test” numbers reported on “2020-07-13”, “2020-07-11”, and “2020-06-12” are not representative of the time trend (reclassification of how tests were reported, dumped on a a few dates) and are omitted
  • the reporting format changed to include “probable” cases on 6/12/2020. A total of 202 cases were classified as “probable” up to that date, and they were all included in the new case count on 6/12/2020.
  • Test reporting Note: Beginning on 6/12/2020 the number of all tests reported includes all tests, including repeat tests in the same person; an unclear number of tests were added on 6/12/2020 to account for this change in reporting.
  • COVID-19 data reported for June 29, 2020 reflect results from 2 days- due to an unplanned shutdown of the state surveillance system during the weekend.
  • Age strata populations are taken from the US Census estimate for TN, year 2019 estimate (the most recent available)
  • Test reporting note: on 4/19/2021 a large number of backlogged results were released and reflected in the 4/17-4/19 results; most cases were not active and reflect results from several weeks before this date.

Vaccinations

A bit of good news. Vaccination data is not updated every day- most recent data provided.

TN now provides detail on number of recipients and number of vaccinations.

Age Group Given 1st Vaccine (n) Estimated Pop. (n) Given 1st Vaccine (% in Age)
<5 10900 408605 2.67%
5-11 119258 585900 20.35%
12-15 156446 343637 45.53%
16-20 230052 428409 53.70%
21-30 548940 943700 58.17%
31-40 573213 879850 65.15%
41-50 589057 854356 68.95%
51-60 677701 893802 75.82%
61-70 682573 782810 87.20%
71-80 452299 481538 93.93%
81+ 191876 226567 84.69%
Population within each age group estimated from US Census estimates for TN, 2019
Vaccination data last updated on 2022-08-22.

most recent County Vaccination data update on 2022-08-22

# max(tndoh_vax_county_current$date)
plot_vacc <-
  tndoh_vax_county_current %>% 
  filter(date == latest_vax_date,
         !is.na(statefp)) %>% 
  ggplot() +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = recipient_count_pct),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f%%", recipient_count_pct), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos), 
                   geometry = geometry, 
                   color = cut(recipient_count_pct, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  labs(fill = "%Pop.:",
       title = paste0("Vaccination Coverage (latest update on ", latest_vax_date, ")"),
       subtitle = sprintf("First Vaccinations: %s (%1.1f%% of TN Population)\nCompleted Vacc. : %s (%1.1f%% of Pop.)\nTotal Vaccinations: %s",
                          comma(max(tndoh_vax_state$recipient_count, na.rm = T), 0),
                          max(tndoh_vax_state$recipient_count_pct, na.rm = T),
                          comma(max(tndoh_vax_state$recip_fully_vacc, na.rm = T), 0),
                          max(tndoh_vax_state$recip_fully_vacc_pct, na.rm = T),
                          comma(max(tndoh_vax_state$vaccine_count, na.rm = T), 0))
  ) +
  theme_void() +
  theme(legend.position = "left") 

# plot_immune <- 
  # tndoh_vax_county_current %>% 
  # filter(date == latest_vax_date,
  #        !is.na(statefp)) %>% 
  # ggplot() +
  # geom_sf(data = TN_counties) +
  # geom_sf(aes(geometry = geometry,
  #             fill = immune_pct),
  #         color = "grey30",
  #         inherit.aes = F) +
  # geom_sf_text(aes(label=sprintf("%0.1f%%", immune_pct), 
  #                  # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos), 
  #                  geometry = geometry, 
  #                  color = cut(immune_pct, 2)), 
  #              size = 3) + 
  # viridis::scale_fill_viridis(option = "C", direction = 1) +
  # scale_color_manual(values = c("white", "black"), guide = F) +
  # labs(fill = "%Pop.:",
  #      title = "Vaccination + Documented Prior Infection") +
  # theme_void() +
  # theme(legend.position = "left") 

plot_vacc

# cowplot::plot_grid(plot_vacc, plot_immune,
#                    nrow = 2)

Rolling 7-Day Averages & Top 10 Most Cases

most recent County data update on 2023-12-09

note: The “New Hospitalizations” in this chart, as reported by the state, may not accurately reflect actual hospitalizations. Right y-axis is population adjusted (per 100k for TN). Ignore for “Test Positivity.”

Counties in the top 10 highest new case rate over the past 7 days; 7-day average trend-line for daily new cases

# add PCT-Pos trend line calculated by STL (Seasonal and Trend decomposition using Loess) 
tndoh_daily$pos_pct_today <- replace_na(tndoh_daily$pos_pct_today, 0)
tndoh_daily$pos_pct_today_tl <- ts(tndoh_daily$pos_pct_today*100, frequency = 7, start = 1)
tndoh_daily$pos_pct_today.stl <- stl(tndoh_daily$pos_pct_today_tl, s.window = 7)$time.series[, "trend"]

# add Vaccine trend line calculated by STL (Seasonal and Trend decomposition using Loess) 
tndoh_daily$new_vaccine_count <- replace_na(tndoh_daily$new_vaccine_count, 0)
tndoh_daily$new_vaccine_count_tl <- ts(tndoh_daily$new_vaccine_count, frequency = 7, start = 1)
tndoh_daily$new_vaccine_count.stl <- stl(tndoh_daily$new_vaccine_count_tl, s.window = 7)$time.series[, "trend"]

variable.labs <- c( "Cases", "Cases", 
                    "Hospitalizations", 
                    "Deaths", "Deaths", "Deaths", 
                    "Tests", "Test Positivity (%)", "Test Positivity (%)", 
                    "Vaccinations", "Vaccinations")
names(variable.labs) <- c( "new_cases_trend.stl", "new_cases_trend.stl", 
                           "new_hosp", 
                           "new_deaths","new_deaths_by_dod", "new_deaths_trend.stl", 
                           "new_tests", "pos_pct_today", "pos_pct_today.stl",
                           "new_vaccine_count", "new_vaccine_count.stl")

plot_state_trends <-
  tndoh_daily %>%  
  mutate(pos_pct_today = pos_pct_today*100) %>%  
  select(date, new_cases_trend.stl, new_deaths_trend.stl, pos_pct_today.stl, new_vaccine_count.stl) %>% 
  pivot_longer(cols = -date,
               names_to = "variable",
               values_to = "value") %>%  
  group_by(variable) %>%  
  arrange(date) %>%  
  ggplot() +
  geom_line(aes(date, value), size = 1, color = "red") +
  labs(x=NULL,
       y="Decomposition 7-day Average",
       title = "TN State Trends") +
  facet_wrap(~variable, scales = "free_y", nrow = 1,
             labeller = labeller(variable = variable.labs)) +
  theme(strip.text = element_text(size = 14))
  
# plot_state_trends <- 
#   tndoh_daily %>%  
#   mutate(pos_pct_today = pos_pct_today*100) %>%  
#   # select(date, new_cases_7davg)
#   # select(date, new_cases, new_deaths_by_dod, new_tests, pos_pct_today) %>%  
#   select(date, new_cases_trend.stl, new_deaths_by_dod, pos_pct_today, new_vaccine_count) %>% 
#   # select(date, new_cases, new_hosp, new_deaths_by_dod, new_tests, pos_pct_today, new_vaccine_count) %>% 
#   pivot_longer(cols = -date,
#                names_to = "variable",
#                values_to = "value") %>%  
#   group_by(variable) %>%  
#   arrange(date) %>%  
#   mutate(value_7davg = RcppRoll::roll_meanr(value, 7)) %>%  
#   ggplot() +
# # plot_state_trends <- 
# #   plot_state_trends +
#   geom_line(aes(date, value), size = 0.6, color = "grey") +
#   geom_line(aes(date, value_7davg), size = 1, color = "red") +
#   # scale_y_continuous(sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="per 100k Population")) +
#   labs(x=NULL,
#        y="Rolling 7-day Average",
#        title = "TN State Trends") +
#   facet_wrap(~variable, scales = "free_y", nrow = 1,
#              labeller = labeller(variable = variable.labs)) +
#   theme(strip.text = element_text(size = 14))


# MIDTN counties plot ####
counties_top10 <-
  tndoh_countynew %>%   #select(new_cases_county_7davg)
  ungroup() %>%  
  filter(date == max(date, na.rm = T),
         name %in% TN_counties$name) %>%  # exclude the pending and out-of-state
  top_n(10, wt = new_cases_county_7davg) %>%  
  pull(name)

plot_top10_counties <-
  tndoh_countynew %>%  
  filter(name %in% counties_top10) %>% 
  ggplot() +
# plot_top10_counties <- 
#   plot_top10_counties +
  # geom_bar(stat = "identity") +
  geom_area(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
  geom_col(aes(date, new_deaths), fill = "red") +
  # geom_line(aes(date, new_cases_county_trend), color = "steelblue", size = 1) +
  geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
  labs(x=NULL,
       y="New Results by day",
       fill = NULL,
       alpha = NULL,
       title = "Top 10 New Cases by Day (7d avg.)") +
  scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
  theme(#axis.text.x =  element_text(angle = 90, hjust = 1),
        legend.position = "none") +
  facet_wrap(~name, nrow = 2, 
             scales = "free_y")+
  theme(strip.text = element_text(size = 14),
        axis.text.x = element_text(angle = 45, hjust = 1))

cowplot::plot_grid(plot_state_trends,
                   plot_top10_counties,
                   nrow = 2, rel_heights = c(1.25,2))

deaths plotted by actual date-of-death. reporting typically lags several weeks, so most recent dates will not be fully reported.

Variants

The Variant distribution is reported by 10 US Regions. Only data for the US total and Region 4 (SouthEast) is shown here. CDC NowCast data is used for this report.

pl_variants <-
  cdc_variants_data %>%
  filter(usa_or_hhsregion == "Region 4" | usa_or_hhsregion == "USA",
         week_ending > lubridate::ymd(Sys.Date())-365,
         published_date == max(published_date)) %>%  
  ggplot(aes(week_ending, share,
             group = variant,
             fill = variant)) +
  ggiraph::geom_area_interactive(color = "grey50", alpha = 0.6, aes(tooltip = variant), hover_css = "stroke:red;stroke-width:inherit;") +
  scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
  scale_x_date(expand = c(0,0), date_breaks = "3 month", date_labels = "%b-%Y" ) +
  theme(legend.position = "none") +
  labs(x = NULL,
       y = "Proportion of cases by Variant",
       title = "Variant Distribution in Southeastern US (CDC Region 4)",
       subtitle = "Hover over with mouse for details",
       fill = "Variant:") +
  facet_wrap(~usa_or_hhsregion)


ggiraph(code = {print(pl_variants)}, zoom_max = 1, width = .8, width_svg = 7, height_svg = 4)

CDC Wastewater Data

Sites in TN recently started reporting data to the National Wastewater Surveillance System (NWSS). (Plots added here Oct 2022).

nwss_dat_tn %>%
  filter(!is.na(pcr_conc_smoothed)) %>% 
  ggplot(aes(date, pcr_conc_smoothed/1000000, group = wwtp_id )) +
  geom_line() +
  facet_wrap(~county_names, nrow = 1, labeller = labeller(county_names = label_wrap_gen(7))) +
  labs(title = "TN Wastewater Surveillance",
       x= NULL, y = "SARS-COV2 Conc by PCR\n(Flow-Population normalized)",
       caption = "Source: CDC, National Wastewater Surveillance System (NWSS)") +
  scale_y_continuous(label=comma, ) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
  theme(strip.text = element_text(size = 14),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position = "none")

Top 10: Largest Increase in Cases

Counties with the top 10 greatest weekly increase in New Cases, by % Change of 7-day average versus 7 days prior (limited to counties with average of >5 daily cases). Case rates are currently increasing in 44 counties.

plot_top10_increase <-
  tndoh_countynew %>%  
  filter(name %in% counties_top10_increase) %>% 
  ggplot() +
# plot_top10_increase + 
  # geom_bar(stat = "identity") +
  geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
  geom_col(aes(date, new_deaths), fill = "red") +
  geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
  labs(x=NULL,
       y="New Results by day",
       fill = NULL,
       alpha = NULL,
       subtitle = "New Cases by Day: 10 Counties with the largest 1-week increase") +
  scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
  theme(#axis.text.x =  element_text(angle = 90, hjust = 1),
        legend.position = "none") +
  facet_wrap(~name, nrow = 2, 
             scales = "free_y")+
  theme(strip.text = element_text(size = 14),
                axis.text.x = element_text(angle = 45, hjust = 1))

plot_top10_increase

# print("No increase in any counties")

Counties at or near all-time peak

# *currently disabled due to drop in cases*  
counties_near_peak_cases <-
  tndoh_countynew %>%   
  ungroup() %>%
  filter(new_cases_county_trend_ratio >= 0.9,
          date >= Sys.Date() - 14,
         !county %in% c("Pending", "Out of State")) %>% 
  group_by(county) %>%
  filter(new_cases_county_trend_ratio == max(new_cases_county_trend_ratio)) %>%
  ungroup()  

counties_near_peak_cases_top10 <-
  counties_near_peak_cases %>% 
  slice_max(order_by =  new_cases_county_trend, n= 10) %>% 
  pull(county) 

No counties are currently near (at or within >90% of) all-time peak new case rate.

Best 10: Largest Decrease in Cases

Counties with the top 10 greatest weekly Decrease in New Cases, by % Change of 7-day average versus 7 days prior (limited to counties with average of >5 daily cases).

Only counties with a decline in case averages are shown (may be <10).
Cases are currently declining in 43 counties

# MIDTN counties plot ####


pl_best10 <-
 tndoh_countynew %>%
  filter(name %in% counties_top10_decrease) %>%
  ggplot() +
  geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
  geom_col(aes(date, new_deaths), fill = "red") +
  geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
  labs(x=NULL,
       y="New Results by day",
       fill = NULL,
       alpha = NULL,
       subtitle = "New Cases by Day: 10 Counties with the largest 1-week Decrease") +
  ylim(0,NA) +
  scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
  theme(#axis.text.x =  element_text(angle = 90, hjust = 1),
        legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1),
        strip.text = element_text(size = 14)) +
  facet_wrap(~name, nrow = 2,
             scales = "free_y")
pl_best10

# cat("There were no counties with a decreased 7d average")

# ifelse(length(counties_top10_decrease) >= 1,
#        # IF >=1 with negative trend
#        plot(pl_best10),
#        # IF No counties with negative trend
#        "There are no counties with a decreasing case trend")

Interactive Tennessee Map

County color indicates the New Cases per day per county population (7 day average)

# make county data to list
county_plots <-
  tndoh_countynew %>% 
  filter(!is.na(date),
         date >= lubridate::ymd("2020-03-15")) %>%  
  group_by(name) %>%  
  nest() %>%  
  arrange(name) %>%  
  mutate(pl_new = map2(data, name,
                       ~ggplot(data = .x) +
                         geom_col(aes(date, new_cases, fill = "Cases"), alpha = 0.5) +
                         geom_col(aes(date, new_deaths, fill = "Deaths")) +
                         geom_line(aes(date, new_cases_county_7davg, fill = "Case 7-Day Avg"), color = "steelblue", size = 1) +
                         # geom_col(aes(date, new_hospitalized, fill = "Hospitalized")) +
                         labs(x=NULL,
                              y="New Results by day",
                              title = glue::glue("{.y} County"),
                              fill = NULL,
                              alpha = NULL,
                              subtitle = "New Results by Date Reported") + 
                         scale_fill_manual(values = c("Cases" = "steelblue",
                                                      "Case 7-Day Avg" = "steelblue",
                                                      "Deaths" = "red"),
                                           # "Hospitalized" = "orange"),
                                           labels = c(paste0("Total Cases: ", comma(.x$total_cases[.x$date == max(.x$date)], digits = 0)),
                                                      paste0("7d Avg: ", comma(.x$new_cases_county_7davg[.x$date == max(.x$date)], digits = 0)),
                                                      paste0("Total Deaths: ", comma(.x$total_deaths[.x$date == max(.x$date)], digits = 0))
                                                      # paste0("Hospitalized: ", .x$total_hospitalized[.x$date == max(.x$date)])
                                           )
                         ) +
                         scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
                         # ylim(0,NA) +
                         scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
                         theme(axis.text.x =  element_text(angle = 45, hjust = 1),
                               legend.position = "bottom")
  ))

county_data_byday_now <-
  left_join(TN_counties, 
            {tndoh_countynew %>%  
                filter(date == max(date, na.rm = T)) %>%   
                select(name, total_cases, new_cases_county_7davg, new_cases_per100kpop_7davg, new_cases_per100kpop_7davg_rank,
                       new_testpos_pct, new_testpos_pct_7davg, new_tests, new_pos_tests,
                       new_cases_rank, new_cases, new_cases_rank, new_cases_county_trend_ratio,
                       total_deaths, total_hospitalized)}) %>%  
  mutate(cases = replace_na(total_cases, 0),
         deaths = replace_na(total_deaths, 0)) %>%  
  left_join({county_plots %>%  select(name, pl_new)},
            by = "name")


lbls <-
  sprintf("<strong><i>%s Co.</i></strong><br/>
          Total Cases: %s<br/>
          New Cases/d/100k, 7dAvg (rank): %1.1f (%s/%s)<br/>
          New Cases Today (rank): %s (%s/%s)<br/>
          Test +%% Today: %1.1f%% (%s/%s)<br/>
          Test +%% (7dAvg): %1.1f%%<br/>
          Total Deaths: %s<br/>
          Current Case trend vs. Max: %1.1f%%<br/>
          <i>Click to view time-course</i>", 
          toupper(county_data_byday_now$name), 
          comma(county_data_byday_now$total_cases, digits = 0), 
          county_data_byday_now$new_cases_per100kpop_7davg, county_data_byday_now$new_cases_per100kpop_7davg_rank, max(county_data_byday_now$new_cases_per100kpop_7davg_rank), 
          comma(county_data_byday_now$new_cases, digits = 0),county_data_byday_now$new_cases_rank, max(county_data_byday_now$new_cases_rank), 
          county_data_byday_now$new_testpos_pct, county_data_byday_now$new_pos_tests, county_data_byday_now$new_tests, 
          county_data_byday_now$new_testpos_pct_7davg,
          county_data_byday_now$total_deaths,
          county_data_byday_now$new_cases_county_trend_ratio*100
          ) %>%  
  lapply(htmltools::HTML)


# pal <- colorQuantile("viridis", domain = county_data_byday_now$new_cases_per100kpop_7davg, n = 10, reverse = FALSE)
pal <- colorNumeric("viridis", domain = county_data_byday_now$new_cases_per100kpop_7davg, reverse = FALSE)


leaflet() %>%  
  addProviderTiles(providers$Stamen.TonerLite) %>% 
  addPolygons(data = TN_counties, # borders of all counties
              color = "grey",
              fill = NA,
              weight = 1) %>% 
  addPolygons(data = county_data_byday_now,
              group = "name",
              fillColor = ~pal(county_data_byday_now$new_cases_per100kpop_7davg),
              fillOpacity = 0.8,
              weight = 1,
              highlight = highlightOptions(
                weight = 5,
                color = "red",
                bringToFront = TRUE),
              label = lbls,
              labelOptions = labelOptions(
                style = list("font-weight" = "normal", padding = "3px 8px"),
                textsize = "15px",
                direction = "auto")) %>% 
  addLegend(pal = pal, 
            values = county_data_byday_now$new_cases_per100kpop_7davg,
            title = "New Cases/100k (7d Avg)",
            opacity = 0.9,
            position = "bottomright") %>%  
  leafpop::addPopupGraphs(county_data_byday_now$pl_new,
                          group = "name",
                          width = 400, height = 200)

TN County Heatmap

Counties are ordered by highest new cases per day (current) based on the most recent 7 day average (per 100k population). Max limit of the color scale set at 200 to increase contrast- some local outbreaks exceed this limit (eg prison outbreaks in Bledsoe, Lake, Wayne Co.)

county_order <-
  tndoh_countynew %>%  
  filter(date == max(date)) %>%  
  mutate(county_bynewcases = reorder(county, new_cases_per100kpop_7davg)) %>%  
  select(county, county_bynewcases)

# no longer interactive due to file size issues.
TN_heatmap_100k <-
  tndoh_countynew %>%  
    filter(!county %in% c("Pending", "Out of State")) %>% 
  left_join(county_order, by = "county") %>%  
  mutate(new_cases = replace_na(new_cases, 0)) %>%  
  ggplot(aes(x = date,
             y = county_bynewcases,
             fill = new_cases_per100kpop_7davg)) +
  geom_tile(aes(fill = new_cases_per100kpop_7davg)) +  
  viridis::scale_fill_viridis(labels = scales::comma,
                                 limits = c(0, n_ul), oob = scales::squish,
                              option = "plasma") +
  viridis::scale_color_viridis(labels = scales::comma,
                               limits = c(0, n_ul), oob = scales::squish,
                               option = "plasma") +
  labs(fill = "New Cases/100k\n(7dAvg.)",
       color = "New Cases/100k\n(7dAvg.)",
       x = NULL, y = NULL) +
  scale_x_date(date_breaks = "2 month", date_labels = "%b-%Y", expand = c(0, 0)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 8),
        panel.grid = element_blank(),
        axis.text.x.bottom = element_text(angle = 90),
        legend.position = "bottom",
        legend.background = element_blank(), 
        panel.border = element_blank(),
        legend.title = element_text(size = 14),
        legend.text =  element_text(size = 10))

# x_TN_heatmap_100k <- girafe(ggobj = TN_heatmap_100k,
#             options = list(
#               opts_hover(css = "stroke-width:0.5")))
TN_heatmap_100k

Notes:

  • 12/20/2020: Increased upper limit to 200 because 125 just wasn’t enough.

New cases by county

Note: assignment of new cases to county usually lags the total count. The State Health Dept. assigns cases to county of residence, rather than testing site. Individual county health departmens may report numbers differently, such as the number of cases diagnosed within the county, regardless of residence

# fig 2 State ####

pl_tnmap_tnsdc <-
  tndoh_countynew %>%  
  filter(!is.na(date)) %>%  
  filter(date == max(date)) %>% 
  mutate(total_cases = case_when(total_cases == 0 ~ NA_real_,
                               TRUE ~ total_cases)) %>%  
  ggplot() +
# pl_tnmap_tnsdc <-
#   pl_tnmap_tnsdc +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = as.numeric(total_cases)),
          inherit.aes = F) +
  geom_sf_text(aes(label= comma(total_cases, digits = 0), 
                   geometry = geometry, 
                   color = cut(total_cases, 2)), 
               size = 2.5) + 
  viridis::scale_fill_viridis(labels = scales::comma,
                                 # limits = c(0, 100), oob = scales::squish,
                              option = "plasma") +
    # viridis::scale_fill_viridis(option = "D", direction = -1) +
    scale_color_manual(values = c( "white", "black"), guide = F) +
  labs(fill = "Total Cases:")+
         # caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via TN State Data Center") +
  theme_void() +
  theme(legend.position = "left")  
# ggsave(here::here("figures/pl_tnmap.jpg"), dpi = 600, width = 12, height = 4, units = "in")

pl_tnmap_7dDelta <-
  tndoh_countynew %>%  #select(new_cases_county_7davg_1wkdelta)
  filter(!is.na(date)) %>%  
  filter(date == max(date)) %>% 
  ggplot() +
# pl_tnmap_7dDelta <-
#   pl_tnmap_7dDelta +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = new_cases_county_7davg_1wkdelta*100),
          inherit.aes = F) +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  geom_sf_text(aes(label=sprintf("%1.1f",new_cases_county_7davg_1wkdelta*100), geometry = geometry),
               # color = cut(new_cases_county_7davg_1wkdelta, 2)),
               size = 3) +
  # viridis::scale_fill_viridis(option = "D", direction = -1) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "1-Week\nChange %:")+
  # caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via TN State Data Center") +
  theme_void() +
  theme(legend.position = "left")  

# pl_tnmap_newhosp <-
#   tndoh_countynew %>%  
#   filter(!is.na(date)) %>%  
#   filter(date == max(date)) %>%  
#   ggplot() +
# # pl_tnmap_newhosp <-
# #   pl_tnmap_newhosp +
#   geom_sf(data = TN_counties) +
#   geom_sf(aes(geometry = geometry,
#               fill = as.numeric(new_hospitalized)),
#           inherit.aes = F) +
#   geom_sf_text(aes(label=new_hospitalized, geometry = geometry), 
#                color = "black",
#                # color = cut(new_hospitalized, 2)), 
#                size = 3) + 
#   scale_fill_gradient(low = "white", high = "red", 
#                       limits = c(0, max(tndoh_countynew$new_hospitalized [tndoh_countynew$date == max(tndoh_countynew$date, na.rm = T)], na.rm = T)))  +
#   scale_color_manual(values = c( "white", "black"), guide = F) +
#   labs(fill = "New\nHospitalizations:") +
#   theme_void() +
#   theme(legend.position = "left")  

pl_tnmap_new_tnsdc <-
  tndoh_countynew %>%  
  # filter(date == max(date)) %>%  
  filter(date == max(date)-1) %>%  
  ggplot() +
# pl_tnmap_new_tnsdc <-
#   pl_tnmap_new_tnsdc +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = as.numeric(new_cases)),
          inherit.aes = F) +
  geom_sf_text(aes(label=new_cases, geometry = geometry), 
               color = "black",
               # color = cut(new_cases, 2)), 
               size = 3) + 
  scale_fill_gradient(low = "white", high = "red", 
                      limits = c(0, max(tndoh_countynew$new_cases[tndoh_countynew$date == max(tndoh_countynew$date, na.rm = T)], na.rm = T)))  +
  labs(fill = "New Cases:") +
       # caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via New York Times") +
  theme_void() +
  theme(legend.position = "left")  
# ggsave(here::here("figures/pl_tnmap_new.jpg"), dpi = 600, width = 12, height = 4, units = "in")

pl_cases_pct <-
  tndoh_countynew %>%    
  # filter(date == max(date, na.rm = T)-1) %>%  
  filter(date == max(date, na.rm = T)) %>%
  ggplot() +
# pl_cases <-
#   pl_cases +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_cases_per100kpop),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f%%", total_cases_per100kpop/1000), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_cases_per100kpop, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Total Affected\n(Cases %Pop.):") +
  theme_void() +
  theme(legend.position = "left") 


pl_FINAL_tnsdc <-
  cowplot::plot_grid(title_tnsdc,
                   pl_tnmap_new_tnsdc,
                   pl_tnmap_7dDelta,
                   # pl_tnmap_newhosp,
                   # pl_tnmap_newtests_tnsdc,
                   pl_tnmap_tnsdc,
                   pl_cases_pct,
                   nrow = 5,
                   rel_heights = c(.1,1,1, 1, 1),
                   align = "y")
  unclassified <- 
    tndoh_countynew %>%  
    filter(date == max(date),
         county_lower %in% c("pending", "out of state"))

  
    # pl_FINAL_tnsdc
  pl_FINAL_tnsdc +
    annotation_custom(gridExtra::tableGrob(
      tibble(
        Category = unclassified$county,
        Cases = unclassified$new_cases), 
      rows=NULL, cols = NULL,
      theme = gridExtra::ttheme_minimal(base_size = 10)),
      xmin=.85, xmax=.95, ymin=0.66, ymax=.9) +
    # annotation_custom(gridExtra::tableGrob(
    #   tibble(
    #     Category = unclassified$county,
    #     Cases = unclassified$new_hospitalized), 
    #   rows=NULL, cols = NULL,
    #   theme = gridExtra::ttheme_minimal(base_size = 10)),
    #   xmin=.85, xmax=.95, ymin=0.31, ymax=.275) +
    annotation_custom(gridExtra::tableGrob(
      tibble(
        Category = unclassified$county,
        Cases = unclassified$total_cases), 
      rows=NULL, cols = NULL,
      theme = gridExtra::ttheme_minimal(base_size = 10)),
      xmin=.85, xmax=.95, ymin=0.26, ymax=.375) 

Deaths and Crude Case Fatality ratio

The estimate of case fatality is presented as the crude deaths/cases ratio, and is not adjusted for age. This could markedly affect the rate in different counties. This is also affected by the ability to test and identify as many cases as possible, which may differ by county.

pl_deaths_total <-
  tndoh_countynew %>%    
  # filter(date == max(date, na.rm = T))  %>%
  filter(date == max(date, na.rm = T)-1)  %>%
  ggplot() +
  # pl_deaths <-
  #   pl_deaths +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_deaths),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.0f", total_deaths), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_deaths, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Total Deaths:") +
  theme_void() +
  theme(legend.position = "left") 

pl_deaths_percapita <-
  tndoh_countynew %>%    
  # filter(date == max(date, na.rm = T))  %>%
  filter(date == max(date, na.rm = T)-1)  %>%
  ggplot() +
  # pl_deaths <-
  #   pl_deaths +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = total_deaths_per100kpop),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.0f", total_deaths_per100kpop), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths_per100kpop, pct_pos), 
                   geometry = geometry, 
                   color = cut(total_deaths_per100kpop, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Total Deaths\nper 100k Pop.:") +
  theme_void() +
  theme(legend.position = "left") 


pl_cfr <-
  tndoh_countynew %>%    
  # filter(date == max(date, na.rm = T)) %>%  
  filter(date == max(date, na.rm = T)-1) %>%
  mutate(cfr = total_deaths/total_cases) %>%  
  ggplot() +
# pl_cfr <-
#   pl_cfr +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = cfr*100),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f %%", cfr*100), 
  # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(cfr, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Case Fatality %:") +
  theme_void() +
  theme(legend.position = "left") 




cowplot::plot_grid(title_tnsdc, 
                   pl_deaths_total,
                   pl_deaths_percapita,
                   pl_cfr,
                   nrow = 4, rel_heights = c(.1,1,1,1))

# ggsave(here::here("figures/pl_deaths_map.jpg"), dpi = 600, width = 12, height = 6, units = "in")

County Population adjusted rates

The estimated Tennessee county populations are taken from the US Census (2019 estimate). The estimate of case fatality is presented as the crude deaths/cases ratio, and is not adjusted for age. This could markedly affect the rate in different counties

pl_newcases_percapita <-
  tndoh_countynew %>%  #select(new_cases_per100kpop_7davg)   
  filter(date == max(date, na.rm = T)) %>%
  # filter(date == max(date, na.rm = T)-1) %>%  
  ggplot() +
# pl_newcases <-
#   pl_newcases +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = new_cases_per100kpop_7davg),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.0f", new_cases_per100kpop_7davg), 
                   # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
                   geometry = geometry, 
                   color = cut(new_cases_per100kpop_7davg, 2)), 
               size = 3) + 
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  # scale_color_manual(values = c( "black", "white"), guide = F) +
  labs(fill = "Daily Cases per\n100K Pop.\n(7d Avg):") +
  theme_void() +
  theme(legend.position = "left") 

# pl_test_percapita <-
#   tndoh_countynew %>%  
#     mutate(total_tests_pctpop = total_tests/population*100000) %>%  
#     # select(total_tests_per100kpop)   
#   filter(date == max(date, na.rm = T)) %>%
#   # filter(date == max(date, na.rm = T)-1) %>%  
#   ggplot() +
#   # pl_cfr <-
#   #   pl_cfr +
#   geom_sf(data = TN_counties) +
#   geom_sf(aes(geometry = geometry,
#               fill = total_tests_pctpop),
#           color = "grey30",
#           inherit.aes = F) +
#   geom_sf_text(aes(label=sprintf("%0.1f ", total_tests_pctpop/100), 
#                    # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos), 
#                    geometry = geometry, 
#                    color = cut(total_tests_pctpop, 2)), 
#                size = 3) + 
#   viridis::scale_fill_viridis(option = "C", direction = 1) +
#   scale_color_manual(values = c("white", "black"), guide = F) +
#   # scale_color_manual(values = c( "black", "white"), guide = F) +
#   labs(fill = "Tests per\n100k Population:") +
#   theme_void() +
#   theme(legend.position = "left") 

pl_test_positivity <-
  tndoh_countynew %>%
  filter(date == max(date, na.rm = T)) %>%
  ggplot() +
  geom_sf(data = TN_counties) +
  geom_sf(aes(geometry = geometry,
              fill = new_testpos_pct_7davg),
          color = "grey30",
          inherit.aes = F) +
  geom_sf_text(aes(label=sprintf("%0.1f ", new_testpos_pct_7davg),
                   geometry = geometry,
                   color = cut(new_testpos_pct_7davg, 2)),
               size = 3) +
  viridis::scale_fill_viridis(option = "C", direction = 1) +
  scale_color_manual(values = c("white", "black"), guide = F) +
  labs(fill = "Test+%\n(7d avg):") +
  theme_void() +
  theme(legend.position = "left")

cowplot::plot_grid(title_tnsdc, pl_newcases_percapita, pl_cases_pct, pl_test_positivity,
                   nrow = 4, rel_heights = c(.1,1,1,1))

# ggsave(here::here("figures/pl_deaths_map.jpg"), dpi = 600, width = 12, height = 6, units = "in")

US Data from the CDC

Case data for these maps are taken from the CDC. and the 2019 state population estimates are taken from the United States Census Bureau.

Formerly I used the COVID Tracking project, but this site is no longer updating data as of March 7, 2021.

CDC Data Currently unavailable

New Cases

mouse over state for details

max_newcaserate <- max(cdc_currentdata_2163$new_cases_state_7davg_100k[cdc_currentdata_2163$date == max(cdc_currentdata_2163$date)], na.rm = T)
                                                  
pl_us_cases_100k <-
  cdc_currentdata_2163 %>%  
  filter(date == max(date)) %>%  
  ggplot(aes(geometry = geometry,
             fill = new_cases_state_7davg_100k)) +
# pl_us_cases_100k <-
#   pl_us_cases_100k +
  ggiraph::geom_sf_interactive(aes(tooltip = tt,
                                   data_id = state), size = .2) +
  viridis::scale_fill_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate), oob = scales::squish,
                              option = "plasma") +
  labs(fill = "Cases/100k/d\n(7dAvg)") +
  theme_void() 

ggiraph::girafe(code = print( pl_us_cases_100k),
                options = list(
                  opts_hover(css = "stroke:orange;stroke-width:2;")))

US Heatmap

max_newcaserate_alldates <- max(cdc_currentdata$new_cases_state_7davg_100k, na.rm = T)

state_order <- 
  cdc_currentdata_2163 %>%  
  filter(date == max(date)) %>%  
  mutate(state_bynewcases = reorder(state, new_cases_state_7davg_100k)) %>%  
  select(state, state_bynewcases)

us_heatmap_100k <-
  cdc_currentdata %>%  
  left_join(state_order, by = "state") %>%  
  filter(date >= lubridate::ymd("2020-04-01")) %>%  
  ggplot(aes(x = date,
             y = state_bynewcases)) +
  geom_tile(aes(fill = new_cases_state_7davg_100k)) +
  viridis::scale_fill_viridis(labels = scales::comma,
                               limits = c(0, max_newcaserate_alldates), oob = scales::squish,
                              option = "plasma") +
  viridis::scale_color_viridis(labels = scales::comma,
                             limits = c(0, max_newcaserate_alldates), oob = scales::squish,
                            option = "plasma") +
  labs(fill = "Cases/100k/d\n(7dAvg)",
       color  = "Cases/100k/d\n(7dAvg)",
       x = NULL, y = NULL) +
  scale_x_date(date_breaks = "1 month", date_labels = "%b", expand = c(0, 0)) +
  theme_minimal() +
  theme(axis.text.y = element_text(size = 12),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.background = element_blank(), 
        panel.border = element_blank(),
        legend.title = element_text(size = 14),
        legend.text =  element_text(size = 12))
us_heatmap_100k

XKCD bonus

You made it to the end, so here’s a random XKCD comic (please excuse any profanity- just randomly selected):

Session Info

Hmisc::markupSpecs$html$session()
 R version 4.3.1 (2023-06-16 ucrt)
 Platform: x86_64-w64-mingw32/x64 (64-bit)
 Running under: Windows 11 x64 (build 22621)
 
 Matrix products: default
 
 
 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base     
 
 other attached packages:
  [1] rvest_1.0.3         httr_1.4.7          Hmisc_5.1-1        
  [4] patchwork_1.1.3     ggtext_0.1.2        htmlTable_2.4.1    
  [7] RcppRoll_0.3.0      cowplot_1.1.1       urbnmapr_0.0.0.9002
 [10] ggiraph_0.8.7       formattable_0.2.1   leafpop_0.1.0      
 [13] leaflet_2.2.0       sf_1.0-14           lubridate_1.9.3    
 [16] forcats_1.0.0       stringr_1.5.0       dplyr_1.1.3        
 [19] purrr_1.0.2         readr_2.1.4         tidyr_1.3.0        
 [22] tibble_3.2.1        ggplot2_3.4.4       tidyverse_2.0.0    
 
To cite R in publications use:

R Core Team (2023). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

To cite the Hmisc package in publications use:

Harrell Jr F (2023). Hmisc: Harrell Miscellaneous. R package version 5.1-1, https://CRAN.R-project.org/package=Hmisc.

To cite the ggplot2 package in publications use:

Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.